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ABSTRACT 


The transient deformations and stresses due to 
earthquakes are studied using time dependent finite element 
models. Efficient computer programs incorporating frontal 
solution and time stepping procedure are developed for the 
modelling of geodynamic problems. This scheme allows us to 
investigate the quasi-static phenomena including the effects 
of the complex rheological structure of a tectonically 
active region. 


From three dimensional models of strike slip 
earthquakes, it is found that lateral variation of viscosity 
affects the characteristics of surface deformations. The 
vertical deformation is especially informative about the 
viscosity structure in a strike slip fault zone. The 
subsiding region near the fault in a layered model can become 
an uplift region with a laterally heterogeneous model. For an 
earthquake with a fault length of 350 km and maximum slip of 
5 m, the post-seismic uplift can reach 15 cm in 25 years if 
the low viscosity zone beneath the fault extends to 20 km deep. 
In such a model, the shear stress in the region adjacent to 
the fault tip along the strike direction increases with time 
rapidly for about two decades following the event, and 
reaches a constant value about 30 years after the event. 

The time dependent stress due to post-seismic relaxation can 
be larger than the stress due to the elastic response. This 
accelerated stress accumulation may have caused the observed 
earthquake migration along the North Anatolian fault zone. 


Fine grid two dimensional models are used to study 
the time dependent deformation following a thrust event in a 
subduction zone. In particular, results from a series of 
models with different viscosity structures are compared with 
the uplift data along a leveling route from Whittier to 
Anchorage following the 1964 Alaskan earthquake. It is ‘found 
that a simple layered model does not explain the data, but 
a model including a descending slab and a low viscosity wedge 
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above the wedge can produce the observed features. 

Three dimensional viscoelastic model of a thrust 
earthquake indicates that the transient disturbance on plate 
velocity due to a great plate boundary earthquake is 
significant at intermediate distances, but becomes barely 
measurable 1000 km away from the source. The accelerated 
stress accumulation along the strike direction following a 
model thrust earthquake is consistent with the observation 
of earthquake migration along the subduction zones. 
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CHAPTER 1 
INTRODUCTION 


In this thesis, we will consider the problem of 
modelling the time dependent deformations and stresses 
associated with earthquakes and aseismic slips. Two and 
three dimensional time dependent finite element schemes have 
been implemented to do the calculations. The purposes of 
this study are to understand the physical processes 
underlying the observed time dependent phenomena of crustal 
movements and seismicity, to develop a scheme of modelling 
that can be used to interpret the data from modern geodetic 
techniques, and to study the rheological structure of the 
earth . 


One of the clearest expressions of current plate 
motions are the great earthquakes that occur at plate 
boundaries. Investigating these events may provide us with 
information about the geodynamic processes in the earth. 
Elastic and viscoelastic motions associated with these 
events are also valuable in probing the rheological 
structure of the earth. Plate tectonics theory has provided 
the basic framework for the interpr etations of earthquake 
mechanisms. Stresses build up as the plates are pulled 
apart, converge together, or slide past each othe*-. 
Earthquakes are simply the fractures to relieve these 
stresses. However, knowledge about many important phases in 
the earthquake mechanisms is still far from complete. 
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In particular, the important problem of the 
accumulation, release, and relaxation of the stresses 
associated with earthquakes is still not fully understood. 
This problem must be considered within the context of the 
tectonics of the region containing the earthquake fault 
zone. A very useful tool for studying this problem is the 
time dependent crustal movements associated with an 
earthquake. These movements result from the interaction of 
the asthenosphere and the lithosphere and provide 
information about the physical processes of stress 
accumulation and relaxation. However, the data are still 
fragmentary; the fault region is usually complex; the 
interpretations of these data are difficult and the 
possibilities of the application to practical problems like 
earthquake prediction have not been fully exploited. 

At the present time, there are intensive surveying 
activities in various seismically active areas in the world. 
There are also efforts to use space technologies for 
geodetic measurements to a precision never before achieved. 
Precise and frequent geodetic measurements covering large 
areas may be available in the near future. The 
interpr etation and modelling of these high quality data will 
require considerable efforts. 

In this thesis we use the finite element method to 
study the time dependent deformations and stress diffusions 
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associated with earthquakes and aseismic slips. An 
earthquake is modelled as an abrupt (step function in time) 
slip and an aseismic slip event is modelled as a prescribed 
time dependent slip. The faulting interacts with the 
viscoelasticity of a region and causes time dependent 
deformations and stresses, and these effects are evaluated 
in the models. When data are available, the model results 
are compared with geodetic and seismicity observations to 
obtain constraints on viscosity structure. The finite 
element scheme is capable of realistically modeling a 
complex region like the subduction zone. Both linear and 
nonlinear rheologies can be used. The models are 
constructed from our knowledge of the regional tectonics, 
and the viscosity distribution in a region is modelled in a 
realistic way. Fine grid two dimensional models are 
constructed to study the effects of vertical and lateral 
variations of viscosity in detail. In particular, two 
dimensional models have been used to study the post-seismic 
uplift of the 1964 Alaskan earthquake. When two dimensional 
models are not adequate to describe the problem, three 
dimensional models are utilized, and the effects particular 
to three dimensional models are evaluated. These include 
the stress diffusions associated with finite strike slip and 
dip slip events, the effects of lateral heterogeneities 
across a strike slip fault zone, and the time dependent 
movements off the symmetry plane of a dip slip event. We ' 
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have adopted the frontal solution technique into the time 
dependent finite element solution, which greatly reduces the 
amount of computer memory needed. The programs can be 
implemented on a moderate size computer while providing 
efficienct calculation speed. 

Chapter 2 discusses the numerical scheme used in 
this thesis. Its merits and shortcomings are compared with 
other solution techniques. The principles underlying the 
frontal solution technique and the extension to time 
depenedent solutions are presented. Test examples are 
given, and the model results using linear and nonlinear flow 
laws are compared . 

In chapter 3 we discuss the time dependent 
deformation and stress relaxation after strike slip 
earthquakes. Three dimensional models with lateral 
heterogeneities are used. It is found that the lateral 
heterogeneities can alter the characteristics of the time 
dependent deformations in a distinctive way. The possible 
relation of time dependent stress and earthquake migration 
along transform fault zone are also discussed. 

Chapter 4 is concerned with the modelling of dip 
slip earthquakes using two dimensional plane strain models. 
The effects of the layered structure, the descending 
lithosphere, and the low viscosities on the crustal 
movements are discussed. In section 4.5, a model of the 
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great 1964 Alaskan earthquake is constructed. The model 
results are compared with the repeated leveling data 
collected after the earthquake. Discussion on the results 
in chapter 4 are given in section 4.6. In particular, the 
various possible mechanisms to explain the observed crustal 
movements after the 1964 Alaskan earthquake are discussed. 
Three dimensional effects of dip slip events are discussed 
in chapter 5, where we consider the perturbation on plate 
motion at far distances by a great dip slip event. Stress 
diffusion phenomenon and its relation to seismicity in 
subduction zones are also discussed. Chapter 6 summarizes 
the principal results and conclusions of this thesis. 
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Chapter 2 

COMPUTATIONAL METHODS FOR TIME DEPENDENT PROBLEMS IN 

GEODYNAMICS 

2.1 Analytical models of crustal movements of earthquakes 

In solving simple problems, analytical methods are 
usually available and desirable, since they give exact 
solutions, are usually simple in form, and the physical 
meaning of each term in the solution can be identified. 
However, in the calculation of time dependent stress and 
deformation in the lithosphere, we soon realize that 
analytical methods are of limited use. Too many 
simplifications must be imposed and a realistic solution is 
usually impossible. Some previous work will be briefly 
reviewed in this section, as they provide a physical picture 
for the process going on in time dependent phenomena, and 
provide check of the validity of numerical solutions. 

The simplest model for an earthquake fault is a 
dislocation in a half space. A strike slip fault is 
represented by a screw dislocation while a dip slip is 
modelled by an edge dislocation. Analytical solutions for 
dislocations in continuum have long been of interest in 
engineering and solid state physics. Steketee (1958a, b) 
first applied general Volterra dislocations to the problem 
of faulting. This formulation has been further pursued by 
Chinnery (1961,1963,1966), Maruyama (1963,1964,1966) and 
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Press (1965). Static dislocation models in a half space 
have been generalized to three dimensions by Mansinha and 
Smylie (1971), and to the spherical earth problem by Ben- 
Menahem et al (1969) and Singh and Ben-Menahem (1969). 

These models have also been generalized to layered media 
(Singh 1970, Sato 1971, Sato and Matsu'ura 1973, Javanovich 
et al 1974ab, Brown 1975), and to general polygonal 
dislocations (Brown 1975). Static dislocations have also 
been extensively used to fit observed geodetic deformation 
by earthquakes (Chinnery 1961, Press 1965, Savage and Hastie 
1966, Hastie and Savage 1970, Fitch and Scholz 1971, Canitez 
and Toksoz 1972, Thatcher 1975ab, Ando 1975, Freund and 
Barnett 1976, Savage 1979). A review of dislocations in 
seismology was given by Savage (1979). 

When repeated survey data are compared with models, 
time dependent effects should be included in the modelling. 
Although the elastic properties of the earth are slowly 
varying functions of space and a simple analytic solution in 
a homogeneous medium may be adequate, the viscosity 
structure of the earth changes by orders of magnitude with 
space. This should be taken into account when modelling the 
time dependent phenomena. 

For a model of a two dimensional homogeneous 
viscoelastic half space, there is no time dependent 
deformation associated with either a screw or an edge 
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dislocation. This can be seen directly from the elastic 
half space solution and the elastic-viscoelastic 
correspondence principle. Briefly stated, the elastic- 
viscoelastic correspondence principle says that if the 
solution of an elastic problem is known, then the Laplace 
transform of the solution to the corresponding viscoelastic 
problem can be found by replacing the elastic constants by 
some quotient of operator polynomials, and the actual loads 
by their Laplace transforms. The time domain solution can 
then be obtained through inverse Laplace transform 
(Christensen 1971). For example, the vertical displacement 
U on free surface for a vertically dipping fault is 


U = 


IT 


2 72 

x + d 


2 2 ' 


( 2 . 1.1 


where U is the fault offset, x is the horizontal 
o 

coordinate, and D and d are the lower and upper limits of 
the fault depth respectively. The surface deformation 
depends only on the geometry of the fault, not on the 
elastic constants. From the elastic-viscoelastic 
correspondence principle, the solutions are thus identical 
for both elastic and viscoelastic solutions. In general, 
surface deformation from a two dimensional dislocation in a 
homogeneous half space does not depend on the elastic 
constants (Savage and Prescott 1978b). Inhomogeneities a^e 
required to introduce time dependent behavior in two 
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dimensional problems. Singh and Rosenman (1974) showed that 
the identitiy of elastic and viscoelastc displacement 
solutions persists fo*' a vertical dip slip fai^lt in a 
homogeneous three dimensional half space. However, three 
dimensional strike slip earthquakes in general have time 
dependent deformation. 

Elsasser (1969) first proposed a one dimensional scheme 
to model the stress diffusion in the lithosphere caused by 
an earthquake. The model is an elastic beam overlying a 
layer of Newtonian fluid, where the elastic beam simulates 
the lithosphere and the viscous fluid simulates the 
asthenosphere . A sudden (step function) push or pull on the 
edge of the elastic plate simulates an earthquake. The 
equation governing the displacement in space and time can be 
reduced to a diffusion equation, and the solution for a step 
function input to such a system is the complementary error 
function. Similar one dimensional models have been further 
discussed by Bott and Dean (1973) and Anderson (1975). 

Melosh (1976) extended this linear one dimensional model 
from Newtonian fluid to power law behavior and applied it to 
the study of the aftershock sequence of the Rat Island 
earthquake of 1965. The common feature of these elastic 
plate over viscous fluid systems is that they have a rather 
prominent stress diffusion effect. The stress at a point 
far away from the source increases rather rapidly with time 
as the stress propagates through the point. Savage and 
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Prescott (1978b) pointed out that this system of elastic 
plate and viscous fluid substratum is not a suitable model 
for studying transient phenomenon. In such a system, 
viscous fluid cannot be deformed instantaneously by shear, 
and the deformation at the time of the earthquake at any 
place other than the source is zero. Because of this, the 
stress caused by the disturbance at any place has a rather 
sharp rise time. This system has been extended to 
viscoelastic media by Ho and Smith (1979). 

Rosenman and Singh (1973ab) obtained the time 
dependent stress and deformation due to three dimensional 
strike slip earthquakes in a homogeneous viscoelastic half 
space by using the correspondence principle and the known 
elastic solution for a strike slip fault. Nur and Mavko 
(1974) first found a solution to large time dependent 
deformation associated with a dislocation in a vertically 
inhomogeneous medium. The model is an elastic layer 
overlying a viscoelastic half space. They used the method 
of images to obtain an approximate solution of a two 
dimensional dislocation in an elastic layer over an elastic 
half space with different elastic constants. The time 
dependent solution of a linear viscoelastic medium can be 
obtained from elastic-viscoelastic correspond ence principle 
Nur and Mavko pointed out that a large thrust type 
earthquake may provide a tool for exploring the rheology of 
the earth's upper mantle. Barker (1976) used Haskell's 
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(1953) method of propagator matrix to obtain a solution of a 
screw dislocation in a layered elastic medium; then he 
obtained the time dependent solution of a viscoelastic 
medium by using the elastic-viscoelastic correspondence 
principle. Rundle and Jackson (1977ab) and Rundle (1978) 
obtained solutions for a strike slip fault in an elastic 
layer overlying a viscoelastic half space. They used the 
method of images to obtain the elastic solution; then they 
used the elastic-viscoelastic correspond ence principle to 
obtain the time dependent solutions. Thatcher and Rundle 
(1979) also obtained a solution for a two dimensional edge 
dislocation in an elastic layer over a viscoelastic half 
space. However, their solution does not agree with that of 
Nur and Mavko. Smith (1974) pointed out that Nur and 
Mavko's solution for an edge dislocation involves some gross 
approximations. Thatcher and Rundle (1979) also pointed out 
that Nur and Mavko’s solution for an edge dislocation cannot 
be reduced to the correspond ing elastic solution at the time 
immediately after the earthquake. Savage and Prescott 
(1978b) extended the Nur and Mavko solution for a screw 
dislocation to the study of the earthquake cycle for strike 
slip faults. Spence and Turcotte (1979) obtained the 
response of an infinite sequence of step-wise offsets on a 
two dimensional strike slip fault. Cohen (1980ab) extended 
Rundle's solution of a three dimensional strike slip 
earthquake to a standard linear solid medium. 
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The analytical solutions for time dependent 
calculations usually use the elastic-viscoelastic 
correspondence principle. The medium property is thus 
restricted to linear viscoelastic materials. At the present 
time, most solutions are restricted to a flat layer over 
half space models. This may be adequate in modelling 
general effects at large distance from a zone of complicated 
tectonics where the earthquake occurs (if three dimensional 
models are available) . But it certainly has limitations in 
applying the layered solution to a complex region like 
subduction zone. Furthermore, the inverse Laplace 
transformation must often be done in approximation, thus the 
exactness of analytical solutions, their other advantage, is 
also lost. 

To model the earth realistically, including 
complications due to structure and materials, we decided to 
take a numerical approach. When the geometry is 
complicated, the numerical approach is not only more 
versatile than the analytic solutions, but is also more 
convenient to use. Once the programs are set up, they can 
be used to compute different solutions; only some input 
parameters have to be changed. A brief discussion on 
numerical methods used in this thesis is given in the 
following sections. 
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2.2 Remarks on available numerical schemes 

There are several well developed numerical methods 
to calculate stresses and strains. The most commonly used 
numerical techniques for solving partial differential 
equations are probably the finite difference method (FDM), 
the finite element method (FEM), and the boundary integral 
equation method (BIE). Historically, the developement of 
each of these three methods has its own origin and physical 
basis. However, as each method developed, the definition of 
each method became more general and the techniques became 
more involved, and the distinction among these three methods 
became less clear. For example, some finite element worker 
interpreted BIE as a kind of finite element process with a 
special form of trial function (Zienkiewicz 1977) * There 
are also finite difference equations derived via finite 
elements (Wempner 1971). If we disregard these artifacts 
due to definitions and consider the basic form of the 
method, each method does have its own merits and 
shortcomings. The choice of the numerical method depends 
very much on the problems to be solved. 


FDM has been used longer than the other two methods. 
The basic essence of this technique is to replace the 
govering differential equations and boundary conditions by a 
set of algebraic equations using finite difference 
approximations. If the boundaries in the problem can be 
described easily with some mathematical formulation, FDM can 


21 


provide a very efficient solution. In such cases the length 
of the program code and computational time are often shorter 
than those of the corresponding FEM solutions. On the other 
hand, it is awkward to use FDM if the boundaries in the 
problem have complex geometric shapes. FDM does not have to 
use regular mesh. points, tut the advantage of simplicity is 
lost when irregular mesh points are used. 

FEM is often (but not always) derived from the 
governing variational principle and Rayleigh-Ritz principle. 
The governing differential equations do not appear in the 
formulation. Each element can be considered as a physical 
element. The behavior of the element depends only on its 
material properties and nodal values. The solution region 
can be built up with finite elements as if the actual object 
is built up by pieces of structural elements. There is 
great freedom in choosing the element type and size to suit 
the actual problem. Complicated geometry of boundaries and 
material inhomogeneities can be handled with ease. In 
particular, the grid space can be designed so that high 
resolution is obtained where it is needed, and coarse grid 
spacing can be used in regions of little interest. On the 
other hand, there are lengthy bookkeeping jobs in finite 
element programs. The computational costs are often higher 
than those of finite difference solutions (if FDM can also 
do the job). Nowadays, both FDM and FEM are standard tools 
in engineering analysis and there are excellent textbooks 
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covering these subjects (e.g. Hildebrand 1968, Zienkiewicz 

1977) , 

The boundary integral solutions are obtained by 
transforming the differential equations governing the 
behavior inside and on the boundary of the domain to 
integral equations over the boundaries. There was no 
general way to solve these integral equations before the 
advent of computers. The use of boundary integrals was, to 
a great extent, limited to theoretical investigations of the 
existence and uniqueness of the solutions. The developement 
of BIE as a real tool for problem solving using computers is 
relatively new, but in recent years it has gained growing 
attention (e.g., Massonet 1965, Mendelson 1973, Lachat and 
Watson 1976, Zienkiewicz et al 1977, Cruse 1977, Cole et al 

1978) . There are a number of advantages for BIE to make a 
scheme competitive with FEM and FDM. The most obvious one 
is that it reduces the dimensionality of the problem by one. 
The numerical solution of the integral equations involves 
the modelling of boundary data rather than volume data. For 
three dimensional problems the numerical analysis is 
performed on the two dimensional boundary surface; thus the 
size of the problem is greatly reduced. The procedure can 
be used to calculate the solution at any interior point, 
unlike FEM or FDM which allows calculation only at nodal 
points. Sometimes it also has advantages in problems with 
infinite region or singularities. However, the matrice 
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involved in BIE are not banded or symmetric. This may well 
offset the reduction of problem size, as efficient 
algorithms exist for solving symmetric and banded matrix 
equations. The solutions are limited to problems which have 
suitable Green's functions, thus its use in non-linear 
material or highly heterogeneous bodies is limited. BIE can 
sometimes be coupled to FEM to combine the advantages of 
both methods (Zienkiewicz et al 1977) 


In this thesis, we will use the finite element method 
exclusively. The main reason for this choice of numerical 
scheme is that we would like to have a versatile method. 
Among the above three methods, FEM takes the most effort to 
program, and the computational cost is the highest. But it 
is also the most versatile. The realistic earth is far from 
a homogeneous body, and the regions where earthquakes occur 
are usully more complicated than other regions. If we are 
interested in modelling the time dependent movements 
associated with earthquakes, we must be prepared to use the 
solution method under various conditions. For an earthquake 
in subduction regions, there are descending slabs dipping at 
various angles and low viscosity wedges above the descending 
slabs. In transform fault regions, the shear heating may 
produce lateral heterogeneities. It is essential that our 
approach can handle these features. We also desire a method 
which can handle a wide varieties of material properties. 

For this reason, we choose the more involved but more 
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versatile finite element method as our numerical method 
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2.3 Finite element solutions of geodynamic problems 

Since FEM is suitable for modelling the 
heterogeneous medium, it has been applied to the study in 
geophysics by many investigators. It has been applied to 
study geological folds using a viscous fluid model 
(Dieterich 1969, Parish et al 1976, Cobbold 1977, Woidt and 
Neugebauer 1980), to the studies of regional stress fields 
and tectonics using either creep or static models (Jungels 
and Frazier 1973, Shimazaki 1974ab, Luo 1974, Smith 1974, 
Bischke 1974, Neugebauer and Breitmeyer 1975, Kusznir and 
Bott 1977, Bird 1978, Neugebauer and Sophon 1978, Melosh and 
Raefsky 1978,1979, Slade et al 1979, Kasahara 1978, Kosloff 
1978, Richardson 1978ab, Richardson and Bergman 1979, Kato 
1979), and to studies of seismic wave propagation (Lysmer 
and Drake 1972, McCowan et al 1977, Smith 1975, Smith and 
Bolt 1976, Archuleta and Frazier 1978, Schlue 1979). In 
addition, there is a large amount of literature using the 
finite element method for studies in soil mechanics and rock 
mechanics. The first study of earthquake deformation using 
the finite element method is probably by Jungels and Frazier 
(1973). They used two dimensional static models to study 
the deformation of the 1972 San Fernando earthquake. Smith 
(1974) extended the analysis to a viscoelastic medium and 
time dependent behaviors using two dimensional plane strain 
models. Two dimensional plane strain models has also been 
used by Bischke (1974) , Melosh and Raefsky (1979) and Slade 
et al (1979). 
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The computational scheme used in this thesis have been 
influenced by the work of Smith (1974). We adopted his 
treatment of the fault model and boundary conditions, 
however, we designed a scheme with wider applicability. 

Smith used the elastic-viscoelastic correspondence principle 
and solved a set of static solutions in Laplace domain using 
finite elements. The time dependence was obtained by using 
numerical Laplace transform (Schapery 1961, Adey and Brebbia 
1973). Smith (1974) was mainly interested in two 
dimensional plane strain models and linear materials. An in 
core finite element solver (FEABL) was used. Later Smith 
(1979) extended his programs to three dimensional 
viscoelastic models. However, some out of core solver must 
be designed to take care of the storage problem for three 
dimensional applications-. We wanted to extend the two 
dimensional analysis to three dimensions, and to cover a 
wider class of material properties and modes of fault 
displacements. To achieve these goals, we developed a 
scheme combining the time stepping approach (Zienkiewicz and 
Cormeau 1974) and frontal solution technique (Irons 1970). 
These will be discussed in the next section. 
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2.4 Frontal solution technique and time dependent 
calculation 

2.4.1 The principles of frontal solution and time stepping 
scheme 

The basic finite element method has been covered in 
many textbooks (e.g. Cook 1974, Zienkiewicz 1977), so it 
will not be repeated here. In finite element models of time 
dependent behaviors associated with earthquakes, we have to 
deal with two practical problems: The first is the large 

size (number of degrees of freedom, storage requirement, a’/d 
computational cost) of the problem, and the second is the 
extension to time dependent calculations. We followed the 
time stepping approach of Zienkiewicz and Cormeau (1974) for 
time dependent calculations, and combined the frontal 
solution technique (Irons 1970) to the time stepping 
approach for an efficient solution scheme and great savings 
in computer memory requirements . 

In the time stepping approach for time-dependent 
calculations, the total strain { e } is divided into three 
parts : 

{el = {e} el + {e} cp + {e}° (2.4.1) 

6 X O 

where { e } is the elastic strain, { e } is the initial 

strain, and { e } ^ is the creep strain. Engineering stress 

and strain are used. We use curly brackets to denote 
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vectors and square brackets to denote matrice. Quite 
generally the constitutive law for creep strain rate and 
stress can be put into the form 

U} cp = [F Ha} (2.4.2) 

where dot 'ndicates time rate of change, IT ] is a symmetric 
matrix (may depend on stress) and (a) is the stress. 

The virtual work principle is: 

) ( 6{e} T {or}dV - ) ( 6{u} T {b}dV 
V V 

- ) ( 6tu} T {t}dS = 0 (2.4.3) 

S 

where lb) is the prescribed body force, {t} is the 
prescribed boundary traction, and {u} is the displacement. 

6 indicates variation and T indicates transpose. 

Integration is over the volume V and traction boundary S 
respectively . 

Let 

{e} = [L ] l u} (2.4.4) 

tu} = [N ] {a} (2.4.5) 

where [L] is the operator relating strain and displacement, 
[N] is the interpolation function, and {a} is the nodal 
displacement. Then equation (2.4.3) becomes 
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S{ a) T [ j[N] T CL] T {a} dV - j[N] T {b} dV 
V V 

- $tN] T {t} dS ] = 0 (2.4.6) 

S 

or 

j[B] T {<r}dV - {F} = 0 (2.4.7) 

V 

where [B] = [L][N], {F} = j[N] 1 {b}dV + (CN] T {b}dS 

V S 

Substituting t cr) = [D]{e} el = [D ] ( { e}- { e} cp - te} °) , equation 
(2.4.7) can be put into the standard form: 

C K 3 { a } = {R } (2.4.8) 

where 

CK ] = J[B] T [D][B]dV 
V 

{R} = {F} + j[B] T [D]{e}°dV + J [B ] T [D ] { e} Cp dV 
V V 

Equation (2.4.8) is the set of linear equations to be solved 
using the frontal solution. {e} cp is obtained in a time 
stepping fashion using the equation {6} cp = Cl* Her}. It can 
be shown that this procedure can be applied to the general 
visco-plastic problems (Zienkiewicz and Cormeau, 1974), 
including plasticity and creep problem as two extreme cases. 
The scheme handles a nonlinear creep material in the same 
way as a linear material. 
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In equation 2.4.8, the vectors {a} and { R } have the 
dimension of the degrees of freedom N in the problem, and 
the stiffness matrix {K } has the dimension of N x N. In 
three dimensional problems, it is very common to have 
thousands of degrees of freedom in a model. The problem of 
computer storage then becomes very serious, so it is 
imperative to fully exploit the banded property of the 
stiffness matrix. We have chosen the frontal solution 
approach (Irons 1970, Orringer 1974) to minimize the in core 
storage requirement. In this approach, the familiar 
Gaussian elimination method is still the central process, 
however the assembly and elimination processes are not 
seperated. We accumulate the contribution to the stiffness 
matrix and consistent nodal force vector element by element 
until one diagonal element K nn of the stiffness matrix has 
received all the contr ibutions from the elements connected 
to it. Then the global force-displacement relation 
associated with K 

nn 

N 

t K nj a j = R „ <2- 1 *-9> 

j=0 

is eliminated. The coefficients in this equation are moved 
to outside storage, and their previous locations in core are 
cleared for storing other coefficients. The remaining 


equations a^e updated according to 
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K nn 


(2.4. 10) 


» K. R 

R* * R t - -H* 


K. 


nn 


(2.4.11) 


Since each degree of freedom accepts contributions only 
from elements connected to it, the in core storage is 
usually greatly reduced (if the elements are numbered 
properly) . At the end of the elimination process, the last 
equation will be in the form of: 


«* » ** » 

Y • * * a - R • * * 

n NN a N ■ N 


(2.4. 12) 


which can be readily solved for the displacement a^ a^ is 
then back substituted into the penultimate equation which 
contains only two degrees of freedom. This back- 
substitution process continues in the order exactly opposite 
to that of the assembly process until all the nodal 
displacements are solved. 


Notice in equation 2.4.8, the stiffness matrix [K] is 
the same for each time step. Only the force vector {R} is 
different. So we do not have to perform assembly and 
elimination for each time step; once the reduced stiffness 
matrix elements have been saved on outside storage, all we 
have to do is to update the force vector {R} and back 
substitute. This constitutes tremendous savings when the 
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size of the problem is large. For the three dimensional 
problems we described in chapter 3 and 5, the c.p.u time for 
each subsequent time step was about 5 to 10 percent of the 
first step (which did the assembly and elimination). 

Because of this huge saving of computation time, we have 
refrained from using implicit time stepping schemes (e.g. 
Hughes and Taylor 1978), which may be more stable .,ut 
requires assembly ‘and elimination for each time step. 

2.4.2 Test examples 

Several examples were used to test the finite element 
programs. First we tested the static deformations of a 
vertical dip slip earthquake. The finite element grid for 
this example is given in figure 2.4.1. The model region is 
1600 km long and 800 km deep. The fault lies on the center 
line of the region and has a constant slip 5 m from surface 
to 75 km deep and tapers off to zero at 80 km deep. We have 
used the treatment of Smith (1974) to model the fault. 

There are double nodes along the fault zone. The fault 
offset is prescribed as relative displacement between the 
double nodes, but the absolute locations of the double nodes 
are allowed to move in response to the stress field. In 
this example, we have assumed that the bottom and sides of 
the region, which are far away from the fault, are rigid. 
Since the earthquake is very far from the external 
boundaries, this assumption has 


little effect on the result. 
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The resulting vertical displacement on the free surface of 
the finite element model is given in figure 2.4.2. The 
analytical solution of the vertical displacement on the free 
surface for a vertically dipping fault with uniform slip in 
an infinite homogeneous half space was given in equation 
2.1.1. The resulting analytical solution of vertical 
displacement on free surface for a vertically dipping fault 
from surface to 80 km deep with a uniform slip 5 m is 
superimposed on the finite element solution in figure 2.4.2. 
The difference between the numerical and analytical 
solutions is insignificant. 

The second test example is a fault dipping 45 degrees, 
with a 5 m offset from the surface to 75 km depth, and 
tapers off to zero at 80 km depth. The finite element grid 
is given in figure 2.4.3. The model region is 2560 km long 
and 1160 km deep. The resulting vertical displacement on 
the free surface of the finite element model is given in 
figure 2.4.4. The analytical solution for vertical 
displacement on the ree surface for a 45 degree dipping 
fault with uniform slip is (Jungels and Frazier 1973) 

U 2(2D-x) 2 + 4x(2D-x) + 6x 2 

U = — [ 5 5 (2.4.13) 

1 1 . 3ff (2D-X; + r- 

2D-x 


+ 4tan - ^ ( 


x 


) + 2 + n ] 


where U is the fault offset, D and d are the lower and 
upper limits of the fault depth, x is the horizontal 
coordinate, and the fault intercepts the free surface at x = 
0. The resulting analytical solution of the vertical 
displacement on the free surface for a 45 degree dipping 
fault from surface to 80 km deep with a uniform slip 5 m is 
superimposed on the finite element solution in figure 2.4.4. 
Again the agreement between the analytical and numerical 
solutions is very good. 

To test the time stepping procedure, we calculated the 

creep of a cylinder under internal pressure. The material 

of the cylinder is elastic in bulk and Maxwellian in 

distortion with bulk modulus 1 unit, shear modulus 1 unit, 

and viscosity 1 unit. The outer boundary of the cylinder is 

rigidly fixed. The internal pressure is 10 units. The 

inside radius of the cylinder is 1 unit and the outside 

radius is 2 units. A schematic diagram of the problem and 

the finite element grid are given in figure 2.4.5. Only a 

section of the cylinder is needed because of the symmetry in 

this problem. The finite element grid has 48 elements and 

130 degrees of freedom. The analytical solution for the 

radial stress S and hoop stress S u as functions of radial 

r h 

distance r and time t are (Flugge 1967): 

Xqa 2 (r 2 - b 2 ) v, 

S = -P [ 1 - 5 -- e" At ] (2.4.14) 

r 2Kor 
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S, = - p [ 1 - 


Xqa 2 (r 2 + b 2 ) Xfc 

- e‘ At ] 


2Kb 


2 2 


with 


(2.4. 15) 


X = 


6Kb^ 

6Kfb 2 + q(3a 2 + b 2 ) 


(/ 


f = 


q = 2(/ 


where P is the pressure, K is the bulk modulus, g is the 
shear modulus, (/ is the viscosity, and a and b are the outer 
and inner radii. 

Figure 2.4.6 gives the finite element result of the 
hoop stress superimposed on the analytical solution, and 
figure 2.4.7 gives the results of radial stress. The very 
good fit between the analytical and numerical results 
confirms that the time dependent finite elemen-t procedure we 
adopted works well . 
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2.5 Discussion of non-linear problems 

The flow property of the earth is vitally relevent 
to many fundamental phenomena in geodynamics. For example, 
the mantle convection and the driving force in plate 
tectonics, the strength and the stress state of the 
lithosphere, and the accumulation of strain in a fault zone 
all depend on the flow properties of rocks. Despite great 
advances in recent years in both experimental and 
observational works, some controversy still exists. 

Inferred from the experimental work (Weertman 1970, Carter 
and Ave’Lallemant 1970, Stocker and Ashby 1973, Kohlstedt 
and Goetze 1974, Weertman and Weertman 1975, Carter 1976, 
Weertman 1978, Goetze 1978, Ashby and Verall 1978), the 
probable creep law under mantle conditions should be a non- 
linear power law creep. However, it is difficult to achieve 
the small strain rate of geological processes in the 
laboratory. Inference on the long term behavior of rocks is 
largely obtained from extrapolation. On the other hand, 
analysis of glacial rebound data suggests that the flow law 
in the mantle may be Newtonian (Walcott 1973, Cathles 1975, 
Peltier and Andrews 1975, Peltier 1976). If the time 
dependent crustal movements associated with earthquakes 
behave distinctively different for linear and non-linear 
rheologies, then we can infer the flow law through the 
modelling of crustal movements with linear and nonlinear 
models. In section 2.4 we mentioned that our finite element 
scheme is suitable for both linear and nonlinear models, 
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however available observation may not be sensitive to the 
form of the creep law (Turcotte et al 1973, Parmentier et al 
1976). We conducted a numerical experiment on the behavior 
of postseismic deformations for linear and nonlinear 
rheologies. Model 2A and model 2B have the same elastic 
response, grid structure, and fault configuration, except 
that model 2A follows power law creep and model 2B follows 
Newtonian creep behavior. The fault model is a uniform slip 
of 3 m from surface to 40 km deep and tapers off to zero at 
80 km deep, with a fault dip of 60 degrees. Both models 
have elastic lithospheres, each 80 km thick. Below this 
depth the creep law for model 2A is assumed to be 
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T = 300 + 1579 x [ 1 - exp( - 7.6x10 _6 D ) ] (2.5.3) 


where D is depth in meters. The relations and the constants 
are taken from Ashby and Verall (1978) for the flow 
properties of dry olivine. All the quantities involved are 
in MKS units. The result of vertical displacement changes 
on the free surface at 0.0, 8.0, 25.8 and 114 years after 
the event is given in figure 2.5.1. 

We would like to see if a linear model can also 

produce such a, result. Model 2B assumes a Newtonian flow 

law. The viscosity is assumed to have a layered structure, 

where the lithosphere is again elastic and 80 km thick. 

From 80 to 180 km deep the viscosity is 1 x 10' J poise; from 

180 to 360 km deep it is 2 x IQ^ 1 poise; from 360 to 600 km 

20 

deep it is 0.8 x 10 poise and below 600 km deep it is 2.2 
21 

x 10 poise. Figure 2.5.2 gives the plot for the vertical 
displacement changes at the same times as in figure 2.5.1. 
The results are very similar to model 2A. The difference is 
probably too subtle to be detected by geodetic data at the 
present time. This numerical experiment demonstrates that 
the crustal deformations are not uniquely determined by the 
form of creep law (although they are in general sensitive to 
the values of viscosity, as we will see in later chapters) . 
As far as the time dependent deformation following an 
earthquake is concerned , linear models contain the essential 
elements of time dependent behavior. In the absence of 
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exact information about the material composition, viscosity 
structure, and flow law in different regions, we choose to 
use linear models to study the relaxation effects of 
viscoelsatic materials in the following chapters. 


Figure Captions - Chapter 2 


Figure 2.4.1 The finite element grid fo- computing the 

deformation of a vertically dipping fault. The model 
fault intercepts the surface at coordinate (0,0). The 
fault extends from 0 to 80 km deep. The grid has 1050 
d.o.f and 496 elements . The model region is 1600 km 
long and 800 km deep. 

Figure 2.4.2 Solution for the vertical displacement of the 
vertical dipping fault. Solid curve is the analytical 
solution from Jungels and Frazier (1973). The symbols 
are the finite element solution using the grid in 
figure 2.4.1. The fault in the analytic solution has 
uniform fault slip 5 m from surface to 80 km deep. The 
fault in the finite element solution has a 5 m fault 
slip from the surface to 75 km deep and tapers off to 
zero at 80 km deep. 

Figure 2.4.3 The finite element grid for computing the 

deformation of a 45 degrees dipping fault. The model 
fault intercepts the surface at coordinate (0,0). The 
fault extends from 0 to 80 km deep. The grid has 1050 
d.o.f and 496 elements. The model region is 2560 km 
long and 1160 km deep. 

Figure 2.4.4 Solution for the vertical displacement of the 
45 degrees dipping fault. Solid curve is the 
analytical solution from Jungels and Frazier (1973). 
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The symbols are the finite element solution from the 
grid in figure 2.4.3. The fault in the analytic 
solution has a uniform fault slip 5 m from surface to 
80 km deep. The fault in the finite element solution 
has a 5 m fault slip from surface to 75 km deep and 
tapers off to zero at 80 km deep. 

Figure 2.4.5 (a) The test problenrof cylinder creep under 

internal pressure. The cylinder is elastic in bulk and 
Maxwellian viscoelastic in distortion. The outer 
boundary is rigidly fixed and the inner boundary is 
under 10 unit of pressure, (b) The finite element 
grid used to calculate the time dependent stress in the 
cylinder creep problem. Because of symmetry, only a 
small section of the cyliner is used in actual 
computation. The finite element grid has 48 elements 
and 130 degrees of freedom. 

Figure 2.4.6 Hoop stress as a function of distance in the 
cylinder creep problem. Curves 1, 2, 3 and 4 are 
analytical solutions at times = 0.0, 1.0, 3.0 and 12 
relaxation times. Symbols are the correspond ing finite 
element solutions. 

Figure 2.4.7 Radial stress of the cylinder creep problem. 
Curves 1, 2, 3 and 4 are analytical solutions at times 
= 0.0, 1.0, 3.0 and 12 relaxation times. Symbols are 
the finite element solutions. 
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Figure 2.5.1 Postseismic vertical displacement change on 
the free surface due to the relaxation for model 2A, 
which assumes a power law creep behavior. The curves 
1, 2, 3, and 4 a r e at time 0.0. 8.0, 25.8 and 114 years 
after the event. 

Figure 2.5.2 Postseismic vertical displacement change on 
the free surface due to the relaxation for model 2B, 
which assumes a Newtonian creep behavior. The curves 
1, 2, 3, and 4 are at time 0.0, 8.0, 25.8 and 114 years 
after the event. 
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1. INTRODUCTION 

Great damage has been caused by shallow strike slip 
earthquakes that occur along plate boundaries in various 
parts of the world. The mechanism of these earthquakes 
has long interested seismologists. The study of geodetic 
measurements of the 1906 San Francisco earthquake led to 
the formulation of elastic rebound theory (Reid, 1910) , 
and this theory has remained a basic tenet of the earthquake 
mechanism. The continuous study of accumulation, release 
and relaxation of stresses near the fault zone has provided 
a more detailed mechanism of strike slip earthquakes (e.g. 
Chinnery, 1961; Scholz and Fitch, 1969; Turcotte and Spence, 
1974; Savage, 1975; Thatcher, 1975a,b; Budiansky and Amazigo, 
1976; Rundle and Jackson, 1977a, b; Savage and Prescott, 1978; 
Savage, 1979; Thatcher, 1979; Turcotte et al. , 1979). Much 
of the study used geodetic measurements near the fault zone. 
In particular, static elasticity and dislocation theory have 
often been used to study the stress and displacement field 
caused by strike slip earthquakes. 

Stress accumulation, release and relaxation are time 
dependent phenomena. This is evident from geodetic data, 
and follows from the migration behavior of earthquakes and 
asthenospheric viscosity. In recent years , there have been 
intensive geodetic and creep measurements in the San Andreas 
fault zone, and ultra precise space technology has been 
applied to geodetic measurements (e.g., Niell et al. , 1979; 
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Smith et^ al_. , 1979). Accurate data will be available in 
the near future on the time dependence of crustal deformation. 
A detailed, three-dimensional time dependent model may be 
necessary for the interpretation of such data. On the other 
hand, earthquake migration phenomenon have been observed along 
plate boundaries, most noticeably along the North Anatolian 
fault (e.g. Mogi, 1968; Allen, 1969; Dewey, 1976; Toksdz 
et al. , 1979) . Explaining this time dependent phenomenon 
also requires time dependent models . 

We calculated the long term time dependent response of 
a set of models of strike slip events. The effect of 
relaxation is isolated in these calculations. Most of the 
previous attempts to model time dependent tectonic phenomonon 
after earthquakes used two-dimensional models (e.g., Nur 
and Mavko, 1974; Bischke, 1974; Smith, 1974; Savage and 
Prescott, 1978; Thatcher and Rundle, 1979; Melosh and Raefsky, 
1979) or simple layer and half space solutions (Rosenman and 
Singh, 1973a, b; Barker, 1976; Rundle and Jackson, 1977a, b,; 
Cohen, 1980a, b; Lehner et al. , 1979). 

Two-dimensional models assume an infinite long fault, and the 
effects in the region beyond the fault tip cannot be described. 
However, in this paper we show that there are significant 
effects in the region beyond the fault tip for strike slip 
events. Laterally homogeneous models that assume no lateral 
heterogeneities across the fault zone oversimplify the 
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near-source problem. Most data indicate the presence of 
lateral heterogeneities near the fault (e.g. Sass and 
Lachenbruch 1973, Zandt, 1978) . 

In this paper, we present time dependent calculations 
for finite strike-slip faults in laterally heterogeneous 
media. We use the three-dimensional finite element method 
to model strike slip events. The forward problem is set up 
to predict time dependent deformation and stress for years 
and tens of years following a modelled earthquake. The 
models are representative earthquakes to show the characteristic 
time dependent features of strike slip events. The boundaries 
of inhomogeneities in the models are kept geometrically simple. 
The model results indicate that geodetic measurements after 
an event may provide information on rheological properties 
near the fault zone which are vitally related to earthquake 
occurrence. 

II. MODEL DESCRIPTIONS 

Fault Models 

We present the computation for two classes of models, 
one for a great earthquake, and the other for a moderate 
size earthquake. Due to uncertainty in the viscosity structure, 
several sets of viscosity values for the same fault model 
are used to show a range of possible results. We are 
interested in the general behavior of relaxation following 
the model earthquake. No attempt is made to model a specific 
region in detail. However, the fault displacement and 
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dimension for great earthquake models (Gl, G2 and G3) are 
comparable to those of the 1906 San Francisco earthquake or 
the 1939 North Anatolian earthquake. The dimension for fault 
models Ml and M2 are appropriate for a magnitude 5. 5-6.0 
earthquake such as the Coyote Lake earthquake of 1979 (M = 5.7) 
(Lee et al. 1979) . 

We model the strike slip earthquake as a sudden slip 
on two sides of a rectangular vertical fault surface. The 
offset values are prescribed , and the subsequent displacements 
and relaxation are computed (details given in the section on 
computational scheme and Chapter 2 ) . We assume that the 
slip on fault stops after the step function slip. This 
probably will happen if deviatoric tectonic stresses are 
relieved near the fault zone and friction again takes over, 
and the region deforms in a coherent manner. However, there 
may be cases where this condition is violated. Slow 
afterslip is sometimes observed after an earthquake (e.g., 
Burford, 1972; Bucknam et al . , 1977; Coppersmith et al. , 

1979) . If afterslip does not last long relative to 
relaxation time (years) the long term effect will be 
similar to a step function. 

The model region for models Gl , G2 and G3 is 2740 km 
long, 2320 km wide and 700 km deep. The fault is a 
rectangular region 350 km long and extends vertically from 
0 to 40 km depth. The relative fault offset is 5 meters 
strike slip, except that near the edges of the fault area 
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it tapers off. The offset tapers off linearly from 5 meters 
at 20 km depth to zero at 40 km depth, and it tapers off 
linearly to zero at 35 km from the tips along the strike 
direction as shown in Figure 1. 

The model region for models Ml and M2 is 196 km long, 

169 km wide and 8Q km deep. The fault is a vertical rectangular 
region 20 km long and 12 km deep. The model earthquake has no 
surface rupture. It has relative offset of 30 cm from 3 to 
9 km depth, and the offset tapers off linearly to zero from 
9 to 12 km depth and 3 to 0 km depth. It also tapers off 
linearly to zero at 2.5 km from tips along the strike 
direction (Figure 2) . 

Material Properties 

There is controversy over the laws governing the creep 
behavior of earth materials (Weertman, 1978) . Linear Newtonian 
behavior of the mantle fits the post-glacial rebound data 
(e.g., Cathles, 1975; Peltier and Andrews, 1976), but laboratory 
rock mechanics experiments exhibit non-linear creep behavior. 

For calculation of perturbation caused by the earthquake, we 
choose the linear viscoelastic model, simply because visco- 
elastic materials contain the essential elements of the time 
dependent relaxation phenomena. Furthermore the assumption 
of linearity greatly simplifies the physical picture, as 
we can separate the effect of perturbations caused by the 
earthquake. We concentrate on the calculation of the relaxation 
effect following strike-slip events; other effects can be 
superimposed due to linearity. 
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The material in the model is assumed to be elastic in 
bulk and Maxwellian viscoelastic in distortion. The short 
term elastic constants of earth vary slowly in space 
(Hadden and Bullen, 1969) . On the other hand, the viscosity 
value changes by orders of magnitude from the lithosphere 
to asthenosphere (Cathles, 1975; Peltier and Andrews, 1976). 
The contribution to viscoelastic relaxation due to changes 
in elastic parameters is therefore relatively unimportant. 

In order that we do not unnecessarily complicate the physical 
picture, we assume all the models have the same instantaneous 
response (elastic) properties of bulk modulus 1.3 x 10 J “ 6 
dyne/cm^ and Poisson's ratio 0.25, which are about the average 
values for crust and upper mantle from seismic studies. 

The different model results will be due to different viscosity 
structures in the models. 

The viscosity of the earth is not a well constrained 

quantity, especially near a tectonically active zone such 

as a transform fault. The lithosphere in general can 

withstand aeviatoric stresses for a million years or longer. 

The thickness of the elastic lithosphere in continental 

regions is probably greater than 50 km. The asthenosphere 

has a viscosity value orders of magnitude smaller than the 

lithosphere. From analyses of post-glacial rebound data, 

the low viscosity layer beneath the lithosphere is probably 

20 

on the order of 10 poise (Cathles, 1975; Peltier and 
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Andrews, 1976). This magnitude of viscosity implies that 
deviatoric stresses cannot be sustained there for a time 
scale larger than several years to decades. 

Near an active transform fault, where only shallow 
earthquakes are observed suggests lateral heterogeneities 
may exist. For example, in the San Andreas fault zone 
earthquakes are usually shallower than 10 to 15 km, 
indicating that stress below this depth is being relieved 
continuously below this depth. This thickness of the 

t 1 ' 

"seismo-genic" layer is simply too small compared to the 
generally accepted lithospheric thickness, suggesting 
the existence of low viscosity zone, plastic yield or 
slow slip at depth below. Lachenbruch and Sass (1973, 

1979) found that the San Andreas fault system is contained 
in a broad zone of high heat flow anomaly. They concluded 
that the thin seismo-genic layer is more brittle than the 
layer beneath it, implying the possibility of a shallow 
low viscosity zone there. Theoretically, the shearing 
motion should also cause temperature and viscosity structure 
to vary laterally (Yuan et al. , 1978) . Three dimensional 
inhomogeneities in elastic parameters are also observed in 
tectonically active regions (Zandt, 1978) although it is 
difficult to estimate viscosity structures from elastic 
parameters . 


There is little information on the viscosity of a 
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fault zone. Budiansky and Amazigo (1976) estimate the 

2 1 

"effective viscosity" of lithosphere to be 10 poise in 
California. Nnr and Mavko (1974) and Smith (1974) analyzed 
the vertical deformation of the 1946 Nankaido earthquake and 
conclude the viscosity value below the elastic lithosphere 
is on the order of lO 1 ^ to 10 20 poise. Thatcher and Rundle 
(1979) found that a viscosity value of about 5 x 10 poise 
in the asthenosphere beneath Japan fits geodetic measurements. 

It is plausible that the shearing motion and higher temperature 

*• »'*■ 

in the fault zone make the viscosity lower than in its adjacent 
regions. We will carry out calculations for a range of 
viscosity values. 

Great earthquakes model G1 is the control model; a 
layered structure is assumed. The lithosphere extends to 
80 km depth; it is given a viscosity value of 10 25 poi.se 
from 0 to 40 km, arid 10 poise from 40 to 80 km depth. 

A low viscosity layer extends from 80 to 180 km depth, 

t , o n 

with a viscosity value of 10 AU poise. Below 180 km, the 

. 22 
mantle viscosity is assumed to be 10 poise as shown m 

Figure 3. 

More realistic models incorporate lateral heterogeneities 
across the fault. If a shallow low viscosity zone exists 
beneath a fault, we may assume it extends to the bottom of 
the seismogenic layer; however, its width is not known. 

If the width is much larger than the fault length, the 
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effects of such a zone will be similar to that of a low 

viscosity layer. We use rather wide low viscosity zones. 

A very narrow low viscosity zone may produce different 

results, probably similar to that of an elastic layer. In 

model G2, a low viscosity zone rising to 20 km below the 

surface and extending 140 km on each side of the fault is 

assumed. The viscosity value is assumed to be 10^0 poise, 

the same as in the low viscosity asthenosphere extending from 

80 to 180 km depth. Calculations show that differences in 

resulting time dependent deformation and stress relaxation 

between model G1 and G2 are significant. Model G3 has 

properties intermediate between model 1 and 2, with the 

low viscosity layer extending to 40 km depth. For moderate 

earthquake models Ml and M2 we assume that under the fault 

the low viscosity zone extends to shallower depths. Far 

away from the fault, the lithosphere is 80 km thick. Both 

models assume that from 12 km to 46 km depth there is a 

low viscosity zone 40 km wide on each side of the fault? 

from 46 km to 80 km depth, it is 70 km wide as shown in 

Figure 4. Fast relaxation time is assumed for model Ml 
19 

(viscosity 10 poise) to establish the maximum possible 
effect of relaxation. Model M2 assumes the low viscosity 
value to be 10^0 poise. 


Computational Scheme 


We use three-dimensional time dependent finite element 
models to calculate time dependent motions following a 
modeled earthquake. We combine the frontal solution technique 
(Irons, 1970) to the unified time stepping approach 
(Zienkiewicz and Cormeau, 1974) for a versatile and efficient 
solution scheme. The calculation scheme is described in 
Chapter 2. 

Due to symmetry of a vertical strike slip earthquake, 
only a quarter of the region is needed to be modelled in 
the numerical scheme. All the models in this study use 720 
elements and 2982 degrees of freedom. The grid structure of 
the moderate earthquake models is a scaled grid of great 
earthquake models. Figure 5 is the top view of the grid 
structure of both classes of models. The three-dimensonal 
model is made up of seven identical grid surfaces. The 
element is in 8 node, 24 degrees of freedom hexahedron 
with 8 gaussian integration stations. Elements with 27 
integration stations have been used on some models. The 
results are nearly identical with those of 8 integration 
stations . 

Since only a quarter of the region is used in the 
numerical scheme, the boundary conditions must simulate 
those of the complete strike slip fault. The fault area 
is contained in the boundary plane x = 0 (Figure 5) . The 
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boundary plane y = 0 is a symmetry surface bisecting the 
physical fault model. On plane y = 0, the displacement 
in x direction is set to zero; on plane x = 0, on the 
fault surface the y displacement is set to half the fault 
offset value and on the rest of the x = o plane the y 
displacements are set to zero. The other external 
boundaries are far away from the fault zone. The effect 
caused by the event decreases with distance from the fault. 

We choose a region large enough such that boundaries far 
away from fault zone have little influence on the behavior 
of the region we are interested in. We used rigid and free 
boundary conditions on the sides far away from the fault zone. 
We present results only for those elements where resulting 
stresses differ by less than a few per cent for these two 
extreme cases. The lower boundary of the region is prescribed 
to be rigid. We used 12 time steps for models Gl, G2 and G3 
to calculate time dependent values for up to 49 years after 
the event. For thin lithosphere and fast relaxing models Ml 
and M2, we used 14 time steps for 9.5 years after the event. 
III. Model Results 

Great Earthquake Models 

The results of time dependent deformation and stress 
from different models are presented and compared in this 
section. Model Gl is the control model with a laterally 
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homogeneous layered structure. The horizontal displacements 
on free surface at selected locations for model G1 are 
shown in Figure 6a. The four figures give the horizontal 
displacements at time = 0, 9.5, 25 and 49 years after event, 
respectively. The initial (time = 0) pattern is typical 
of a strike slip fault in elastic media (cf. Chinnery, 

1961) . Afterwards the displacement increases with time 
although the displacement rate decreases . The initial 
displacements decrease very fast with distance from the 
fault. The relaxation process spreads the deformation 
outward from the fault zone. 

In contrast with model Gl, model G2 has a low 
viscosity zone extending to shallow depth near fault. 

The horizontal displacements on the free surface at selected 
locations for model G2 are given in Figure 6b. The 
instantaneous response is the same as model Gl, since they 
have the same elastic parameters. However, the magnitude 
of time dependent displacement is much larger and 
concentrated near the fault zone (where the low viscosity 
zone is shallow) . The time dependent effects can be seen 
more clearly in Figures 7a and 7b, where the displacements 
along the line perpendicular to the center of fault 
(x axis in Figure 1) are shown. Near the fault zone the 
time dependent deformation is in general small compared to 
the instantaneous response for model Gl, while the shallow 
low viscosity zone in model G2 significantly increases 
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the magnitude of time dependent displacement. 

The contours of vertical displacement on free surface 
immediately after the event are shown in Figure 8a. Also 
shown (Fig. 8b) are the shear stress contours. The results 
are what we expect from a shallow strike slip event. In 
the vertical displacement figure there is subsidence in a 
broad area in the upper right quadrant, except near the 
fault tip. Our grid resolution cannot resolve the very 
small uplift zone near the fault tip, but the overall 

r - 

pattern is very similar to that of the analytic half space 
solution of Chinnery (1961) . Subsequent time dependences 
show characteristic differences between layered and lateral 
inhomogeneous models. The contours of vertical displacement 
change after the event (total displacement minus 
instantaneous response) at 5 and 25 years after the event 
for model G1 and are given in Figures 9a and 9b. The 
results of model G2 are given in Figures 10a and 10b. On 
the upper right quadrant, for model G1 the time dependent 
vertical movement is continuous subsidence near the fault 
zone, and uplift away from it. This result is expected 
provided that the horizontal displacement is spreading out 
from the near fault region. The magnitude of this vertical 
movement is on the order of several centimeters. In 
contrast, the results for model G2 is continuous uplift 
throughout the upper right quadrant. The physical reason 
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for this difference can be seen from the horizontal 
displacement plots (Figures 6a and 6b) . In the upper right 
quadrant in the figure, the "flow pattern" of the strike 
slip event is such that material enters the near fault zone 
parallel to the strike direction, but leaves the fault zone 
in the direction 45 degrees from the strike direction at 
large distances. The low viscosity channel near the fault 
zone makes it easier for material to enter the fault zone; 

the thicker lithosphere beyond the low viscosity channel 
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forms a barrier for material to fan out, thus material 
piles up near the fault zone and a bulge results. The 
uplift in model G2 is about 15 cm in 25 years. This effect 
is measureable by geodetic means , and could serve as a tool 
for investigating the viscosity structure of the fault zone 
Another interesting phenomenon is the time dependent 
character of the stress relaxation for these two different 
models. The instantaneous perturbation of horizontal shear 
stress component a xy at 10 km depth caused by the strike 
slip event is shown in Figure 8b, and a X y at 25 years later 
is shown in Figures 11a and lib for models G1 and G2 
respectively. Before an event, due to relative plate 
motion the stress component c? x y should be positive near the 
fault zone for a left lateral strike slip fault as in the 
models. The earthquake relieves the prestress along the 
fault, but reinforces the prestress in front of the fault 
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tip. More detailed time dependency can be seen in Figures 
12a and 12b, where cr^ vs. time at selected positions is 
shown. It can be seen that along the side of the fault, 
cfjjy is negative, implying the prestress is relieved. However, 
viscoelastic relaxation in general accelerates the stress 
recovery. In front of the fault tip, cr xy is positive, and 
the prestress is reinforced. This result of reinforcement 
of prestress has been considered to be the cause of secondary 

faulting (Chinnery, 1966) or creep and aftershocks in front 

/* 

of fault tip (Scholz et al. , 1969) . The earthquake may in 
time trigger subsequent events along the same fault. 

However, for the uniform thick lithosphere model Gl, the 
time dependent stress is small compared to the instantaneous 
response, in this case subsequent earthquakes should happen 
immediately after the event rather than being delayed for 
years. Aftershocks located at the end of the main shock 
fit this picture, however, the short time delay of the 
aftershock is most likely due to local inhomogeneities, 
time dependent friction and creep type relaxation rather 
than large scale mantle relaxation (Dieterich, 1972; 

Mikumo and Miyatake, 1979, Yamashita, 1979). On the other 
hand, for the laterally inhomogeneous model G2, a significant 
portion of the perturbing stress in front of the fault tip 
is accumulated years after the event. The perturbing stress 
levels off a few decades after the event. This accelerated 
stress accumulation in front of the region of a strike slip 
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fault makes the chance of earthquake happening greater during 
the several decades after the event. Triggered events can 
then happen years after the triggering event, as long as 
the stress diffusion is reinforcing the prestress. The 
earthquake sequence after the 1939 North Anatolian event 
(Toksflz et al. , 1979) is consistent with this stress 
diffusion mechanism. It is a general result that the 
relaxation effect near surface becomes more important as 

the thickness of the lithosphere becomes thinner (Rundle 

** 

and Jackson 1977a, b; Savage and Prescott 1978; Cohen, 1980a, b; Lehner 
et al. , 1979). Rundle and Jackson (1977b) and Lehner et al. (1979) 
also pointed out that stress in front of the fault tip 
increases with time after a strike slip event for layered 
models. 

It has been reported (Thatcher 1975a) that the strain 
accumulation accelerates near the fault after the 1906 
San Francisco event, the average strain rate parallel to the 
fault was about 2.5 x 10"® per year for about 25 years after 
the event, and 0.6 x 10"® per year afterwards. For model 
Gl, the strain at a point 18 km away from the fault accumulates 
at an average rate of about 0.1 x 10"® per year, much less 
than the observed value. This simple model does not explain 
the 1906 data. For the laterally inhomogeneous model G2, 
the average strain rate at 18 km away from the fault is 
about 1.0 x 10 ^ per year, which is of the same order of 
magnitude as the .observed values of the 1906 San Francisco 
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earthquake. The model is not a model of the San Francisco 
earthquake in detail, however, the results indicate that if 
the low vi^cosiry extends to shallow depth under the fault, 
viscoelastic relaxation may contribute significantly to the 
time dependent deformation. 

Model G3 has the low viscosity channel to 40 km depth 
instead of 20 km i,'i model G2. The results are intermediate 
between G1 and G2. Figure 13 compares the contours of 
vertical displacement change at 25 years after the event 
for model G2 and G3. The result for model G3 in the upper 
right quadrant is a broad region of uplift, and a small 
region of slight subsidence near the fault tip. This 
indicates that time dependent vertical displacement is 
sensitive to the viscosity structure. Figure 14 compares 
the shear stress o^y at selected positions vs. time for 
models G2 and G3. There is accelerated strain rate a few 
decades after the event for model G3, but the magnitude is 
smaller than that of model G2. 

Moderate Size Strike Slip Earthquake Models 

Models Mi and M 2 are models of moderate size earthquakes 
with buried faults. Model Ml is the fast relaxation model, 
with viscosity 10 ^ poise in the low viscosity zone, while 
model M 2 uses viscosity 10^0 poise instead. The horizontal 
displacements at selected locations on free surface are 
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shown in Figures 15a and 15b. The results are generally a 
gradual increase in magnitude of the displacements. Figure 
16 show’s the instantaneous response of vertical deformation 
on free surface. In the upper right quadrant, there is a 
relatively broad region of uplift near the fault, and 
subsidence far away. This is expected for a fault with 
large fault depth to length ratio (0.6 in this case) . 

The vertical deformation changes at 9.5 years later are 
given in Figures 17a and 17b. For both models Ml and M2, 
in the upper right quadrant, the vertical deformation is 
subsidence near the fault zone and uplift far away from the 
fault, similar to that of the layered medium case. This is 
probably because the width of the low viscosity channel is 
relatively large compared to the earthquake fault dimension. 
The magnitude of time dependent vertical deformation is 
small, generally less than 1 cm. 

Figure 18 gives the time dependence of shear stress 
cr x y at selected locations at 7.5 km depth. We again see 
the stress is relieved along the fault, and reinforced 
in front of the fault. The time dependence is quite 
different for the two models: significant changes are 
completed within 5 years after the event for model Ml, 
while the relaxation is linear in model M2 for the first 
decades. Though the magnitudes of quantities involved in 
those moderate earthquakes are of marginal use with the . 
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precision of present day available geodetic data, in the 
future when more precise and frequent geodetic measurements 
are available, the models may be used to study the 
detailed structure near a fault zone. 

IV. Discussion 

In this study, we calculated the time-dependent stress 
relaxation and deformation for large and moderate size 
earthquakes using different models. To compare these to 
observations we look at space-time migration of seismicity 
and to geodetic data. Migration of seismicity along plate 
boundaries has been observed in South America (Kelleher, 
1972) , the West Pacific (Mogi, 1968) , Alaska-Aleutian 
(Kelleher, 1970) , the San Andreas Fault zone (e.g. Wood 
and Allen, 1973; Lee et al ■ , 1979), and the North Anatolian 
fault zone (e.g. Richter, 1958; Mogi, 1968; Savage, 1971; 
Dewey, 1976; Toksflz et al . , 1979) . In the case of the 
North Anatolian faults there was a relatively quiescent 
period in seismicity prior to 1939. A bi-directional 
trend of seismic migration followed the great earthquake 
of 1939 (Toksflz et al. , 1979) . Savage (1971) proposed a 
dislocation wave theory to explain the earthquake 
migration phenomenon. The occurrence of a sequence of 
events in rapid succession can also be due to an accelerated 
stress accumulation process. In all the model results in 
section III, there is an accelerated stress accumulation 
process in the region in front of the fault tip. The 
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magnitude of this accelerated stress accumulation is 
significant if a low viscosity zone und£r the fault extends 
to a depth of 40 km or less below the surface. It is 
important to clarify that other processes (such as change 
of material properties due to stress level change, slip 
at depth) may contribute to accelerated strain accumulation 
and earthquake migration. These have not been investigated 
in detail. 

A possible scenario for the occurrence of a sequence 
of earthquake follows. The initial earthquake happens when 
stress accumulation exceeds the strength. of the fault zone. 
After this event, the stress accumulation accelerates in 
the region in front of the fault tip. The next earthquake 
happens when the combination of diffused stress and initial 
stress exceeds the strength; in turn this "triggered" 
earthquake triggers the next event in an adjacent region. 
The process continues until the stress along the whole 
fault zone is relieved. This episode is then followed by 
a slow stress accumulation stage and relative quiescence 
in seismicity. Aftershocks that immediately follow the 
earthquake are probably due to local stress adjustments. 
Creep resulting from stress changes at the immediate 
vicinity of the source may result in relatively rapid 
stress adjustments in the source area. Migration of 
earthquakes in time may be related to viscoelastic 
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relaxation and stress diffusion. However, it is not 
possible at this time to test this scenario against other 
possibilities. 

Time-dependent horizontal and vertical motions after a 
strike slip event strongly depend on the viscosities under 
and near the fault zone. Although horizontal displacements 
are much larger than vertical displacement, in this study 
we found that the tiros dependent behavior of the vertical 
displacement is very sensitive to lateral heterogeneities 
of viscosity distribution. The bulge or subsidence formed 
after a great earthquake is of measureable magnitude. Thus 
levelling, in addition to horizontal geodetic measurements, 
after a great strike-slip earthquake may reveal the structure 
near the fault zone. 

In all the models in this study, we have concentrated 
on the relaxation effects after a strike slip event. Some 
possible effects during an earthquake cycle are not explicitly 
modeled. The slow strain accumulation by tectonic stress and 
possible deep slip are not included in the finite element 
calculation. However, since linear material is assumed, 
the effects can be superimposed upon the relaxation process. 

If the recurrence time of the earthquake is on the order of 
hundred years and relaxation time of the asthenosphere is on 
the order of several years, then relaxation after the 
earthquake is a relatively fast process compared with 
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strain accumulation. If we take the state just before the 
earthquake as the reference state, the deformation and stress 
pattern should not be much influenced by the slow strain 
accumulation from the tectonic process . 

Deep aseismic slip below the seismo-genic layer is 
suspected of playing an important role in earthquake 
mechanism (e.g. Thatcher, 19 7 5a, b; Rundle and Jackson, 1977b; 
Savage and Prescott, 1978, Turcotte et al. , 1979). Several 
studies have found that it is difficult to distinguish the 
effects of deep aseismic slip and viscoelastic relaxation 
from geodetic measurements alone (Barker, 1976; Rundle and 
Jackson, 1977a, b; Savage and Prescott, 1978) . This argument 
was based on calculations using laterally homogeneous layered 
models and comparisons of the surface deformation caused by 
viscoelastic relaxation and a given dislocation at depth. 

As pointed out by Savage and Prescott (1978) , there are 
inherent difficulties in using a simple layered model to 
compare these results, because in such cases a distribution 
of slip at depth can always duplicate the viscoelastic 
result. This inherent difficulty may not arise if the 
symmetry is broken by three dimensional inhomogeneity; in 
such cases the images dislocation may not lie on the fault 
surface. The relaxation is more intensive near a low 
viscosity zone, as we have shown in section III. Although 
present day available geodetic data are not constraining, 
geodetic measurements in the future could resolve this question 
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In conclusion, we have implemented a versatile scheme 
to model the time dependent behavior after earthquakes. 
Non-uniform fault slip and three dimensional heterogeneities 
can be included in this scheme. The model results predict 
a stress diffusion phenomenon in front of fault tip after a 
strike slip event: if low viscosity extends to shallow depth 
near the fault zone, the shear stress in front of the fault 
tip will increase significantly with time. The time dependent 
deformation on free surface is more concentrated near the 
fault zone in that case than it is in the case of a laterally 
homogeneous layered structure. The time dependent behavior 
of vertical displacement near the fault may be completely 
altered by the presence of lateral inhomogeneities. 
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Figure Captions 

Fig. 1. Fault for models Gl, G2, and G3. 

(a) Schematic diagrams of fault, double hatched 
area has maximum fault slip, single hatched area 
has tapered fault slip. 

(b) Fault slip along strike direction (Y direction) . 

(c) Fault slip vs. depth (Z direction) . 

Fig. 2. Fault for models Ml and M2 

(a) Schematic diagram of fault, double hatched area has 
maximum fault slip, single hatched area has tapered 
fault slip. 

(b) Fault slip along strike direction (Y direction) . 

(c) Fault slip vs. depth (Z direction) . 

Fig. 3. Sectional views of viscosity distribution for models (a) Gl 
(b) G2 f (c) G3. Numbers with exponent are viscosity in poise. 
Fig. 4. Sectional view of viscosity distribution for model 

Ml. Model M2 has the same structure except that the 

. . 20 . . . 
viscosity is 10 poise m low viscosity zone. 

Fig. 5. Top view of the finite element grid. The three 

dimensional model is made of seven identical plane grids. 

Fig. 6. (a) Horizontal displacements on free surface due to 

strike slip event in model Gl at 0, 9.5, 25 and 49 

years after the event. The location of the fault is 

indicated by a thick line segment and sense of motion 

is indicated by a pair of arrows. 
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Dots are the locations where displacements are 
calculated. Displacement is indicated by a line 
segment from the dot. A scale for the displacement 
(100 cm) and a scale for the map (400 km) are also 
shown in the figure. 

(b) The same plot for model G2. 

Fig. 7. (a) Horizontal displacement vs. time along the line 

perpendicular to the center of fault (x axis) on free 
surface for model Gl. The distance from the center of 
the fault for points 1, 2, 3 , 4, 5, 6 and 7 are 35, 70, 
105, 140, 210 and 280 km. Schematic diagram of the 
locations are shown at the bottom of the graph. 

(b) The same plot for model G2. 

Fig. 8. (a) Contours of vertical displacement on free surface 

immediately after event for great earthquake models 
Gl, G2, and G3 . Broken lines (negative values) indicate 
subsidence; solid lines (positive values) indicate uplift. 
The location of fault is indicated by a thick line segment 
and a pair of arrows. The numbers near the contours are 
uplift or subsidence in mm. The tick marks on the 
frame are at half fault length interval (175 km) . The 
elastic response is the same for model Gl, G2, and G3. 

(b) Contours of shear stress component a X y at 10 km 
depth immediately after the event for models Gl, G2, and 
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G3. The numbers near the contours are stress in bars. 
Stress concentrations are placed near fault tips from 
interpolation, the actual grid resolution is 1/10 of 
the fault length. 

Fig. 9. (a) The vertical displacement on free surface at 

5 years after the event minus the elastic response for 
model Gl. Numbers are amount of uplift or subsidence 
in mm. Ticks are at half fault length intervals. 

(b) The same plot at 25 years after the event. 

Fig. 10. (a) Vertical displacement on free surface at 5 

years after the event minus the elastic response for 
model G2. 

(b) The same plot at 25 years after the event. 

Fig. 11. (a) Contours of shear stress component a X y at 10 km 

depth 25 years after the event for model Gl. 

(b) Same plot for model G2. 

Fig. 12. (a) Shear stress component o X y vs. time at selected 

locations for model Gl. The locations are at 10 km depth. 
A schematic diagram indicating the horizontal positions 
relative to the fault is given to the right of the 
figure. The distances from the fault center along the 
strike direction for points 1, 2, 3, 4, 5, 6, 7, 8, and 
9 are 18, 53, 88, 123, 193, 228, 263, 307 and 385 km. 

The distance from the fault perpendicular to the strike 
direction is 35 km for point 9, 18 km for the rest of 
the points. 

(b) Same plot for model G2. 
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Fig. 13. (a) Contour of vertical displacement on free surface 

at 25 years after the event minus the elastic response 
for model G2 (Identical to figure 10b, included for 
comparison) . 

(b) The same plot for model G3. 

Fig. 14. (a) Shear stress component cr X y vs. time at 

selected locations for model G2. The locations are 
at 10 km depth. A schematic diagram indicating the 

horizontal positions relative to the fault is given to 

*■ 

the right of the graph (identical to figure 12b, included 
for comparison) . 

(b) The same plot for model G3. 

Fig. 15. (a) Horizontal displacements on free surface due 

to strike slip event in model Ml at 0, 1, 3 and 9.5 years 
after the event. 

(b) The same plot for model M2. 

Fig. 16. Vertical deformation on free surface immediately 
after the event for models Ml and M2. Ticks are at 
half fault length interval (10 km) . Numbers near the 
contours are uplift (positive values) or subsidence 
(negative values) in mm. 

Fig. 17. (a) The vertical displacement on free surface at 

9.5 years after the event minus the elastic response 
for model Ml. Numbers are uplift or subsidence in 
mm. Ticks are at half fault length interval (10 km) . 
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(b) The same plot for model M2. 

Fig. 18. (a) Shear stress component o X y vs. time at selected 

locations for model Ml. The locations are at 7.5 km 
depth. A schematic diagram indicating the horizontal 
positions relative to the fault is given to the right 
of the figure. The distance from the fault center along 
the strike direction for points 1, 2, 3, 4, 5, 6, 7, 8, 9 
and 10 is 1.3, 3.8, 6.3, 11.3, 13.8, 16.3, 18.8, 21.9, 

27.5 and 33.8 km. The distance from the fault perpendicular 
to the strike direction is 3.8 km for point 9 and 10, and 
1.3 km for the rest of 'the points. 

(b) The same plot for model M2. 
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Chapter 4 

Two dimensional analysis of time dependent movements 
associated with dip slip earthquakes. 

4 . 1 Introduction 

In this chapter, we will study the time dependent 
deformations and stresses associated with dip slip 
earthquakes using two dimensional finite elements. Various 
fault zone structures and rheological properties are used in 
the calculations. Time dependent stresses and deformations 
for both abrupt and slow aseismic slips are modelled. An 
earthquake is modelled as a sudden slip (step function in 
time) on a dipping fault plane. An aseismic slip event is 
modelled as a prescribed time dependent slip. Viscoelastic 
relaxation is included in the calculation. The medium 
property is assumed to be elastic in bulk and viscoelastic 
in distortion. In this modelling, we are primarily- 
interested in the large earthquakes in subductior, regions. 
These earthquakes result from plate convergence; therefore 
they are very relevant to the dynamics of plate movements. 
Comparison of the model results nf time dependent 
deformations and stresses with geophysical observations may 
help us understand the mechanisms of plate interactions in 
subduction regions. 

A complete description of a region on earth should 
be three dimensional. However, three dimensional finite 
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element models are expensive to implement. The computing 
cost and human efforts in preparing three dimensional models 
can be an order of magnitude more than those of two 
dimensional models. In many cases, two dimensional models 
are adequate to describe the problem. In subduction 
regions, a great earthquake can have a fault length of about 
a thousand km. For example, the 1964 Alaskan earthquake had 
a fault length of about 600 to 800 km (press and Jackson 
1965, Plafker 1969); the I960 Chilean earthquake had a fault 
length of about 1000 km (Plafker 1972). The fault edges 
should have little effect on the surface deformations and 
stresses in the middle part of the fault, which can then be 
calculated using two dimensional models. 

To study the phenomena associated with dip slip 
events, we'll use two dimensional and three dimensional 
models in a complementary way. In a subduction zone, where 
a two dimensional plane strain model is applicable, the 
stresses and deformations of a great dip slip earthquake are 
rapidly varying functions of the distance perpendicular to 
the trench; therefore a fine grid structure is needed to 
delineate these variations. We'll use two dimensional 
models with fine grids to study the phenomena near the fault 
zone. On the other hand, at a distance far away from the 
fault zone, two dimensional models cease to be valid. 
However, the stresses and deformations of an earthquake 
change slowly with space and time at far distances. Thus 




101 


coarse grid three dimensional models can be used. In 
chapter 5, we will use a th^ee dimensional model with 
relatively large grid space to investigate the time 
dependent deformations and stress diffusion phenomena at far 
distances . 

In section 2 of this chapter, we'll briefly review 
the relevant observations in subduction regions. In section 
3 we will discuss some basic results of the relaxation of a 
dip slip earthquake in a viscoelastic medium. In section 4 
we will discuss the interaction of deep aseismic slip and an 
earthquake along different segments of the fault. In 
section 5 we will apply relaxation models to the study of 
post-seismic uplift following the 1964 Alaska earthquake. 

In section 6 we will discuss and summarize the results of 
this chapter. 




4.2 B^ief review of deformation data in subduction regions 


4.2(1 General structure of a subduction region 

Subduction regions are among the most active and 
complex tectonic provinces on earth. They have been 
intensively studied by geophysicists, especially since the 
advent of plate tectonics theory. A comprehensive review is 
beyond the scope of this thesis. We will concentrate on the 
observations which are relevent to time dependent 
deformations. Topographically, a subduction region features 
a deep sea trench, an outer trench bulge, cordillera or 
island arc and sometimes a shallow marginal basin. There 
are volcanoes, heat flow and gravity anomalies. There is 
also intense seismic activity. 85% of the seismic energy 
released in the world is released in subduction regions. It 
is impressive that these complex phenomena are common to 
many subduction zones. This suggests that there are common 
physical processes operating in subduction regions. In 
plate tectonics theory, a subduction region is where a cold 
lithosphere descends into the mantle. All the features of 
subduction regions are probably related to the lithospheric 
subduction process. 

The largest earthquakes in a subduction zone occur 
at shallow depth. To model a g*"eat shallow earthquake, 
attention must be given to the observed features the^e. An 
anomalously low attenuation and high seismic velocity body 


in the mantle is identified to be the sinking lithosphere 
(Isacks et al 1968). The deep sea trench marks the place 
where subduction takes place. With few exceptions (Kanamori 
1971a), great shallow earthquakes are of thrust type. These 
great shallow earthquakes occur along the zone of contact 
between the descending and overiding lithospheres. The 
descending lithosphere acts as a stress guide; intermediate 
and deep focus earthquakes occur inside it (Isacks and 
Molnar 1971, Smith and Toksoz 1972, Toksoz et al 1973). 
Volcanic activity is often present above the sinking 
lithosphere. High heat flow anomalies are o.ften found in 
volcanic arcs and marginal basins while low heat flow 
anomalies are found near deep sea trenches (Uyeda and 
Vacquier 1968, Watanabe et al 1977). The high heat flow in 
volcanic arcs and marginal basins suggests that the 
temperature in the wedge area above a descending slab may be 
relatively high. This is also supported by the observations 
that the attenuations of seismic waves are often high in 
this wedge area (Oliver and Isacks 1967, Utsu 1971, 

Barazangi et al 1975), and that seismicity is absent beyond 
an aseismic front in the wedge (Yoshii 1979). The 
temperature distribution in subduction regions has been 
subjected to intensive investigations (e.g. McKenzie 1969, 
Turcotte and Oxburgh 1969, Minear and Toksoz 1970, Hasebe et 
al 1970, Griggs 1972, Toksoz et al 1973, Andrew and Sleep 
1974, Schubert et al 1975, Toksoz and Bird 1977, Toksoz and 
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Hsui 1978, Hsui and Toksoz 1979, Sydora et al 1980). All 
the evidence indicates that large temperature variations 
exist in subduction regions. There is a cold descending 
lithosphere but a hot wedge area above the slab. It is not 
likely that a simple layered model can adequately describe 
these features in a subdnction region. 

The strain accumulation and relaxation of an 
earthquake are influenced by the viscosity distribution. 

From laboratory study of rock flow properties, the viscosity 
depends strongly on temperature and to a certain extent on 
pressure (Carter 1976, Goetze. 1978, Ashby and Verall 1978). 

A complex thermal regime implies a complex viscosity 
structure. Since viscosity controls the time dependent 
behavior, we will include the viscosity variation of a 
subduction region in our model calculations of time 
dependent phenomena. 


4.2.2. Crustal movements associated with earthquakes in 
subduction regions 


Anomalous crustal deformations before and after an 
earthquake have often been reported (e.g. Imamura 1928ab, 
Tsuboi 1930, Okada and Nagata 1953, Tsubokawa et al 1964, 
Thatcher 1975ab, Fujii and Nakane 1979, Chang 1980). The 
deformations before an earthquake are indicators of strain 
accumulation, and they can serve as premonitary warnings for 
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a coming earthquake. The post-seismic deformations are 
related to the relaxation and redistribution of stresses 
near the fault zone. In this section, we will briefly 
review some observed time dependent crustal movements 
associated with earthquakes in subduction regions. 

Although geodetic data associated with earthquakes 
have been reported for nearly a century, the information is 
still fragmentary, and the phenomena are diverse. Geodetic 
measurements are also prone to errors and at times are 
controversial (e.g. Savage 1975b, Jackson et al 1980). 

There is also the problem of uniqueness in the 
interpretation of geodetic data. At the present time, it is 
not likely that we can interprete the crustal movements in a 
region uniquely. Very often several mechanisms can produce 
the same observed phenomenon. However, geodetic 
measurements do put some constraints on possible models. In 
the future, with the frequent measurements covering broad 
areas using modern techniques like VLBI and laser ranging, 
it may be possible to gain a more thorough understanding of 
the dynamics of plate boundary movements through the 
modelling of geodetic measurements. At the present time, 
the best geodetic data associated with earthquakes in 
subduction zones are either from Japan or from the 1964 
Alaskan earthquake. A review of the geodetic observations 
from these two regions are given in the following 


subsections . 


Information f’-om Japan 


a . 


The largest body of geodetic mearsu**ements 

associated with dip slip earthquakes come f»"om Japan. The 

Geographical Survey Institute of Japan and its predecessors 

carried out first order triangulations and levelings several 

times in the entire country since the last decade of last 

century (Harada 1976). The Geographical Survey Institute of 

Japan made a commendable effort to publish the survey data 

in Journals (Journal of Geodetic Society of Japan, Bulletin 

of Geographical Survey Institute of Japan). However, much 

of the data are still in notebooks and not easily accesible 

to outside users. This one hundred years covered several 

great earthquakes in Japan. As far as time dependent 

deformations are concerned, the most discussed event is 

probably the magnitude 8.2 Nankaido earthquake in 1946 

(Nagata and Okada 1947, Okada 1950, Okada and Nagata 1953, 

Miyata 1955, Fitch and Scholz 1971, Kanamori 1972a, 1973, Nur 

and Mavko 1974, Smith 1974, Bischke 1974, Ando 1975, 

Thatcher and Rundle 1979). This earthquake and the 1944 

Tonankai earthquake (magnitude 8.0) broke much of the 

northern plate boundary between the Philippine plate and the 

Eurasian plate. The source mechanisms for both events were 

determined to be low angle thrust from geodetic and seismic 

studies (Fitch and Scholz 1971, Kanamori 1972a, Ando 1975). 

The combined fault arna of these two events evaluated from 

2 

seismological study was about 10,000 km and the fault slip 
was about 3 meters (Kanamori 1972a) , but the fault area 
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determined from geodetic data was about 25,000 km and fault 
slip was from 3 to 18 m (Fitch and Scholz 1971). 

Kanamori( 1 973 ) interpreted the difference by long lasting 
aseismic and anelastic deformations. Precise levelings were 
carried out several times before and after the earthquake 
(Miyabe 1955). A typical set of geodetic data are shown in 
figures 4.2.1 and 4.2.2. Figure 4.2.1 shows a leveling 
route passing Muroto promontory on Shikoku Island, roughly 
perpendicular to the trench. Figure 4.2.2 shows the 
preseismic, coseismic and postseismic vertical deformations 
along the route (taken from Kanamori 1973). Before the 
earthquake, the trench side subsided with respect to the 
inland. The earthquake reversed the trend by uplifting 
Muroto promontory and subsiding the inland. After the 
event, the post-seismic deformation was uplift near the zero 
crossing line of the coseismic vertical deformation. Fitch 
and Scholz (1971) interpreted the preseismic subsidence near 
Muroto promontory as the overiding plate being dragged down 
by the descending oceanic plate. Thatcher and Rundle (1979) 
proposed that viscoelastic relaxation caused the subsidence. 
For the post seismic deformation, Nur and Mavko (1974) and 
Smith (1974) proposed that viscoelastic relaxation was the 
reason for the deformation. Thatcher and Rundle (1979) 
suggested that deep aseismic reverse slip caused the uplift. 
Fitch and Scholz (1971) proposed that the upper part of the 
fault was undergoing reverse slip while the lower part 
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slipped forward. This really demonstrates the nonunique 
nature of the probelm. At the present time, it is not 
possible to give a unique interpretation . Part of the 
difficulty is that the data were only available on land. 

All the models were forced to fit the data with the trailing 
part of the deformation on land while the largest 
deformation occurred in the sea. To understand the mechanism 
of trench earthquakes, we should seek more independent 
evidence . 

Another often discussed example is the 1923 Kanto 
earthquake. This earthquake was located at the western end 
of Sagami Trough fault, Sagami Trough fault is a transform 
fault between the Asian plate and the Philippine sea plate. 
The levelings in Boso penisula, which started only 20 km 
away from Sagami trough, provided the deformations near the 
source. The fault determined from geodetic and seismic 
studies both indicated that the earthquake fault movement 
was an oblique slip on a dipping fault zone (Kanamori 1971b, 
Ando 197*0. The earthquake fault parameters determined from 
geodetic measurements were 6 m of right lateral strike slip 
and 3 m of reverse dip slip. The dip angle of the fault 
plane is 45 degrees. Thus it is inappropr iate to model this 
fault with a two dimensional dip slip model, especially near 
the end of the fault. The lack of symmetry of the fault 
makes such a three dimensional model expensive. In 
addition, the regional structures are complicated by 
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subsidiary faults and by blocks moving more or less 
independently from each other (Scholz and Kato 1978), making 
the modelling mor# complex. This can be seen from the data 
scattering around the theoretical values from simple models 
(Scholz and Kato 1978, Thatcher and Rundle 1979). 

Two other often discussed earthquakes in Japan as 
far as geodetic measurements are concerned are the 1964 
Niigata earthquake and the 1973 Nemuro-oki earthquake. 

There were premonitary crustal movements before the Niigata 
earthquake detected from geodetic measurements (Tsubokawa et 
al 1964). The fault parameters of the 1964 Niigata 
earthquake have been determined from seismic study by Aki 
(1966) and Abe(1975), and from geodetic study by Abe (1975). 
The fault parameters determined by Abe (1975) were a reverse 
fault dipping 56 degrees with 3.3 m average slip. The fault 
dimension was 80 km long and 30 km wide. Fujii (1978) 
examined the leveling data near Niigata since the turn of 
the century. He found that the precursory anomalous 
deformation near fault zone started about a decade before 
the event. The anomalous crustal movements started from the 
northern part and developed southwards. Near the middle 
section of the fault in the Nezugaseki region on the 
footwall side of the earthquake, the general trend of the 
vertical motion was continuous uplift since the turn of the 
century. The region started an accelerated uplift about one 
decade before the event, became subsidence a few years 
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before the earthquake. Fujii (1978) proposed a creep 
dislocation model to interpret this preseismic anomaly: due 
to tectonic stress, the deep part of the lithosphere (below 
about 50 km depth) started slipping aseismically about one 
decade before the event. This caused a broad uplift zone. 

As the dislocation extended near the free surface, part of 
the uplift zone located on the footwall side began to 
subside. The model in general seems to fit the preseismic 
observation. However, Niigata was a region under rapid 
subsidence because of heavy ground water pumping (Hayashi 
1970), and some levelling results might not represent the 
tectonic movements. 'Even though the recent geodetic data 
can be explained, the significance of this earthquake in the 
context of long term tectonic movements is not certain. 

This earthquake was an intraplate event. It had a rather 
steep dipping fault as compared to the shallow dipping fault 
in Nankaido. Nezugaseki was an area of uplift since late 
Quaternary, as deduced from the heights of coastal terraces 
(Matsuda 1976). The earthquake caused Nezugaseki to 
subside; the subsequent motion of the tidal gage in 
Nezugaseki was continuous subsidence for several years after 
the event (Yohiko Crustal Movement Observatory 1973). 
Nakamura et a.l ( 1 964) suggested that another earthquake 
would uplift the Nezugaseki region to match the coast 
terraces in the future. This earthquake has not happened 
yet . 


Ill 


The Nemuro-oki earthquake of 1973 was a predicted 
earthquake. Five great earthquakes ruptured much of the 
region from Kurile Islands to Hokkaido during the 17 years 
from 1952 to 1969; the only seismic gap in that region 
before 1973 was off the coast of Nemuro plane in eastern 
Hokkaido (Utsu 1972, 1974). The crustal deformations in 
eastern Hokkaido before 1973 were horizontal contractions in 
the direction perpend icular to the trench, and extensive 
subsidence with the subsidence rate increasing from the 
inland area toward the ocean. These were similar to the 
observations before the 1946 Nankaido earthquake, and it was 
consistent with the concept that the overiding plate was 
dragged down by the underthrusting plate (Shimazaki 1974a). 
From the levelling data collected since the turn of the 
century to 1973, the maximum subsidence rate was about 10 
mm/year in Eastern Hokkaido (Tada 1974, Shimazaki 1974ab, 

Abe 1977a). However, from the fossil shell bed along the 
coast, the inferred subsidence is about 2 meters during the 
past five thousand years (Geological Survey of Japan 1974). 
Also from the height change of the marine terraces for the 
past half million years, the subsidence rate is only about 
0.4mm/year (Shimazaki 1974a), an order of magnitude smaller 
than the subsidence rate in this century. If the sea level 
has remained more or less constant (Fairbridge 1961, Curray 
et al 1970), the rapid subsidence in Eastern Hokkaido for 
the past 70 years and the seismic gap strongly suggested a 
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coming earthquake (Utsu, 1972, Shimazaki 1974ab) . The 
predicted earthquake occurred on June 17, 1973. The fault 
mechanism and location were exactly as predicted. The fault 
was determined to be a low angle thrust fault from both 
seismic and geodetic data (Shimazaki 1974b, Tada 1974). The 
aftershock zone almost exactly filled the gap left by the 
five previous great earthquakes (Shimazaki 1974b) . But the 
magnitude of the earthquake (7.4) was smaller than expected. 
Based on the long term subsidence rate, the coast of eastern 
Hokkaido was expected to rebound upward following the 
earthquake (Kasahara 1975). However, the tidal stations at 
Hanasaki and Kushi*-o, which are located on the coast of 
eastern Hokkaido, did not rebound upward after the event 
(Abe 1977). Thus the significance of this earthquake is not 
certain yet. 

The cases of Nemuro-oki and Niigata earthquake 
demonstrated some curious facts about the relationship 
between recent geodetic history and the long term behavior 
in Japan. In regions of recent subduction like southwestern 
Japan, where there are no deep earthquakes and no developed 
trenches, the seismic deformations conform to the Quaternary 
vertical deformation. Consistency of geological and 
geodetic deformations are also found in Kanto and Nobi 
basins, where rapid subsidence is happening (Matsuda 1976). 
On the other hand, in Northeastern Honshu and Hokkaido, 
where there are deep seismicity, developed trenches, 
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marginal basins and a long history of plate subduction, the 
recent geodetic deformation is not consistent with the 
Quaternary geological deformation. Kato (1979) compiled a 
map for the vertical deformation for the past seventy years 
in northern Honshu; in general there is a broad region of 
subsidence, and the subsidence is largest near the coast 
(about 30 cm in 75 years). However, northern Honshu was an 
uplift region during the Quaternary (Matsuda 1976). The 
area of subsidence is too extensive to be treated as a 
premonitary phenomenon for an earthquake in a specific 
region. The great Sanriku earthquake of 1933 happened off 
shore from this region, and it did not reverse the trend of 
the subsidence. Maybe in these regions, a great earthquake 
does not necessarily signal the beginning of a cycle of 
crustal movements, 70 years of geodetic measurements might 
not reveal all the phases of the crustal motion there. The 
genesis of the structures in these well developed trench- 
arc-back arc regions is not fully resolved yet (Uyeda 1977, 
Toksoz and Bird 1977, Hsui and Toksoz 1979, Toksoz and Hsui 
1978). The forces acting in these regions may be complex, 
and this may well influence the geodetic observations. 

To avoid the structural complexity and possible 
contaminations of geodetic data by artificial causes like 
industrial ground water pumping, we will not model the 
geodetic observations in Japan. Instead, we will look into 
the data from the great 1964 Alaskan earthquake, which is 
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reviewed next. 

b. The 1964 Alaskan earthquake 

Geodetic measurements in the subduction regions 
outside Japan usually do not have such a long history. 
However, two great subuction zone earthquakes which happened 
in the last two decades were outside Japan. They were the 
magnitude 8.5 Chilean earthquake of I960, and the magnitude 
8.2 Alaskan earthquake of 1964. Of these two events, there 
are more geodetic measurements associated with the 1964 
Alaskan earthquake. Repeated surveys have been made since 
1964 and may provide useful information (Brown et al 1977). 
Since this earthquake happened in the recent past, we will 
construct a model and calculate the time dependent 
deformation associated with it (section 4.5). It is hoped 
that future surveys will continue and new data can then be 
compared with the model results. 

The 1964 Alaskan earthquake happened in the 
northeastern end of the Aleutian arc, where the arc 
gradually merges eastward into a zone of shallow seismicity 
and transform faulting (Plafker 1969). Surface deformation 
caused by this event has been reported by various 
investigators (Meade 1969, Parkin 1969, Plafker 1969, Small 
and Wharton 1969). The aftershock region in the continental 
shelf and the Gulf of Alaska has been surveyed by echo 
sounding (Huene et al 1967). The regional crustal movements 
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of the 1964 Alaskan earthquake are more extensive than any 
known historic seismic events. The significant tectonic 
deformations, involving uplift, subsidence and horizontal 

p 

displacements, affected a minimum area of 200,000 km . The 
earthquake and its aftershocks were located between the 
trench and the volcanic arc in Alaska. All the larger 
aftershocks (magnitude >5.0) were shallower than 40 km 
(Algermissen 1965). Subsidiary faulting also occurred 
within the deformed region related to the Alaskan 
earthquake. The vertical deformation along a representative 
profile perpendicular to the strike of the trench is given 
in figure 4.2.3 (from Hastie and Savage 1970). There was a 
broad zone of hinterland subsidence, and the maximum 
subsidence was 2.3 meters. Very large uplift was found on 
Montague Island. The zones of uplift and subsidence were 
seperated by a line of zero land level change without abrupt 
offset (figure 4.2.3). The maximum uplift was 11.3 m on 
Montague Island. However, the uplift on Montague Island was 
probably associated with the Patton Bay fault, which is a 
subsidiary fault about 35 km long on land and extends to the 
Gulf of Alaska. Geological relations across the fault 
suggested that this fault was not a major tectonic boundary 
(Plafker 1969). The closest leveling data to the trench was 
on Middleton Island about 60 km off the trench; the vertical 
displacement there was about 4 m. The fault dip was 
determined to be 5 to 15 degrees from a body wave fault 
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plane solution (Stauder and Bollinger 1968), about 20 
degrees from a surface wave nodal plane solution (Kanamori 
1970), and only 3.7 degrees from geodetic data (Hastie and 
Savage 1970). Figure 4.2.3 also shows the theoretical 
surface deformation in the optimal dislocation model by 
Hastie and Savage (1970). The fault parameters given by 
Hastie and Savage had a reverse dip slip of 12.2 m, and a 
left lateral strike slip coimponent of 9.7 m. The fault 
length was 600 km and width 204 km. The dip angle was 3.7 
degrees and the depth of the top of the fault was 4.8 km. A 
secondary fault near Montague Island was included in the 
model to match the abrupt uplift on Montague Island. This 
secondary fault was 300 km long and the dip angle 37 degrees 
with a reverse slip component 17 m and left lateral strike 
slip 12 m. The secondary fault was considered to make an 
important contribution to the observed displacements only 
within 40 km of the fault trace (Hastie and Savage 1970). 

As can be seen from figure 4.2.3, the fit of observations 
with the dislocation model is not perfect. However, as 
pointed by Hastie and Savage (1970), a more complex model 
probably will add too much complexity without much physical 
significance. There are few data near the trench. It is 
not certain if the earthquake fault reached the surface 


there . 
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The coseismic deformation and geological record on 
the coast in Alaska in general were consistent with the 
concept of lithosphere subduction in plate tectonic theory. 
The Alaskan earthquake is the megathrust between the 
descending and overiding lithospheres. On Middleton Island, 
there were steplike flights of marine terraces. A 4 m 
terrace was identified as the marine surface uplifted during 
the 1964 earthquake (Plafker 1969). Radiocarbon dating for 
the older terraces suggests that the last previous uplift 
was about 1400 years ago, and the intervals for sudden 
uplifting was about 500 to 1400 years for the last several 
thousand years (Plafker 1969). The crustal movements along 
shorelines, were generally consistent with the coseismic 
movement. Much of the zone which uplifted and part of the 
zone which subsided during the earthquake showed evidences 
of submergence for the past thousand years. This was 
interpreted as the subsidence caused by the dragging of the 
descending lithosphere before the earthquake (Plafker 1972), 
as in the case of the 1946 Nankaido earthquake. 

The National Geodetic Survey performed first order 
leveling surveys from Whittier to Anchorage in 1964, 1965, 
1968 and 1975. The results were given by Brown et al 
(1977). Figure 4.2.4 shows a simplified map of south 
central Alaska and the leveling route, and figure 4.2.5 
shows the vertical deformation after the 1964 earthquake. 

The route from Whittier to Anchorage was approximately 
perpendicular to the trench strike. Figure 4.2.5 indicates 
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that a bulge was formed along the route. The leveling was 
tied to the tidal gage in Anchorage. The precision was 
believed to be better than 1 . 0 mm per km (Brown et al 1977). 
The uncertainty accumulates proportional to the square root 
of distance. The maximum uplift at the 1975 leveling 
occurred at bench mark D73, which is 40 km to Anchorage; the 
uplift was 54.8 cm there and the standard deviation was 
claimed to be less than 1 cm there. In the 1975 leveling, 
the route was extended from Anchorage to Palmer along a 
route roughly parallel to the trench (figure 4.2.4). There 
was no bulge along this route. This suggested the two 
dimensional nature of the deforamtion. From the studies in 
chapter 3 on the time dependent movements after the strike 
slip earthquakes, we know that the vertical deformation 
caused by strike slip motion is small near the middle part 
of the fault. We can thus use two dimensional plane strain 
models to study the time dependent uplift along this route. 
These data provided a rare set of precision leveling 
covering a long route. In the section 4.5, we will 
construct a model of the 1964 Alaskan earthquake and compare 
the model time dependent behavior with this data set. 

Before we construct a model for the 1964 Alaskan 
earthquake in section 4.5, we'd like to look into some basic 
physical processes of the relaxation after an earthquake 
using simple models. This is given in section 4.3, where we 
discuss the effects of viscosity on the deformation pattern 


after an earthquake. 
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4.3 Effects of layered structure on postseismic deformation 

In this section we look into some basic features of 
the postseismic deformation of two dimensional dip slip 
earthquake models. Smith (1974) studied extensively the 
relaxation behaviors after thrust earthquakes. Our study in 
this section closely follows his approach. Time dependent 
deforamtions after two dimensional dip slip events have also 
been studied by Nur and Mavko (1974) and Thatcher and Rundle 
(1979) using layer over half space models. The purpose of 
this section is to demonstrate the effects of increasing 
viscoity with depth on relaxation, and to discuss the basic 
physical processes in viscoelastic interac tions . 

Smith (1974) found that vertical deformation 
following a dip slip earthquake is informative about the 
fault zone structure. He found that the ratio of fault 
depth to lithosphere thickness is an important controlling 
factor of the behavior of post seismic deformation. If the 
fault does not fracture the whole lithosphere, the area near 
the fault tip will subside after the earthquake. The 
results in this section are consistent with his findings. 
However, there are other parameters which will also affect 
the mode of relaxation. In order to understand the physical 
processes in relaxation, we use only simple models in this 
section. Throughout this chapter we will assume that the 
material is elastic in bulk and Maxwellian viscoelastic in 
distortion. We assume that the medium has a simple 
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instantaneous response of Poisson's ratio 0.25 and Young's 

12 2 

modulus 2.0 x 10 dyne/cm , the same values used in chapter 

3. The first model (model 4A) has a buried shallow dipping 

fault. The schematic diagram of the fault and the finite 

element mesh used in this model is given in figure 4.3*1. 

In the coordinate system of the finite element grid, the 

fault is located between the horizontal coordinates -200 km 

and 0 km, and between 5 and 40 km in depth. The fault is a 

thrust fault dipping 10 degrees. The maximum fault slip is 

assumed to be 12 m, and it tapers off to zero along the 

broken line near the ends in figure 4.3.1. The viscosity 

structure of model 4A simulates a flat elastic lithosphere 

overlying a viscoelastic half space. The lithosphere is 80 

25 

km thick and has viscosity 10 poise. For transient 

phenomenon lasting for decades the response of this high 

viscosity material is practically elastic. The 

20 

viscoeleastic half space has a viscosity 10 poise. We 
calculate the displacement vectors at selected location;] at 
time = 0.0, 1.1, 3*6 and 11 years after the earthquake. The 
snapshots of the deformation are given in figures 4. 3* 2a to 
4. 3. 2d. In these figures, the small circles indicate the 
locations where the calculations are made. The line segment 
stemming from a circle indicates the displacement vector at 
that location. The scale for the displacement is given in 
the lower left corner while the scale for the map is given 
in the lower right corner. The thick line indicates the 
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location of the fault. Figure 4.3.2a is the elastic 
deformation associated with the earthquake. Due to the 
presence of the free surface, the deformation pattern 
differs from the symmetric pattern of a dislocation in full 
space; however, the basic deformation pattern still retains 
the features of a double couple source. The displacements 
above the fault zone in general are larger than those below 

the fault. On the free surface, there is a broad zone of 

> 

subsidence near the lower part of the fault. The subsidence 
changes continuously to uplift toward the upper part of the 
fault. The maximum uplift occurs near the upper tip of the 
fault. Beyond the upper tip of the fault, the vertical 
surface deformation suddenly drops to nsar zero. This 
feature is well known for shallow dipping thrust faults 
(Mansinha and Smylie 1971). 

Although the largest coseismic deformation and 
stresses are located near the fault, the most important 
cause of postseismic deformations lies in the low viscosity 
region below the lithosphere. The elastic lithosphere 
cannot have time dependent deformation by itself. When the 
low viscosity region starts to flow in response to stresses, 
the elastic lithosphere will deform to maintain continuity 
and equilibrium. The region below the lower tip of the 
fault is under compression after the earthquake, as can be 
seen easily from the deformation pattern. Thus material 
tends to move away from that region. In the regions on 
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either side of the fault, the material at far distances is 
drawn towards the fault zone because of the tension created 
by the earthquake. Thus the basic time dependent pattern of 
material movement is that material moves away from the 
compressed region below the fault, and moves towards the 
fault zone from the two sides. This phenomenon can be seen 
more clearly from the snapshots of the postseismic 
displacement changes (total displacement minus coseismic 
displacement) in figures 4.3.3a to 4.3.3c, which show the 
postseismic displacement changes at 1.1, 3.6 and 11 years at 
the selected locations. Two circuits of material movements 
are formed after the earthquake. The resulting vertical 
deformation changes on the f’-ee surface are shown in figure 
4.3.4. There is a subsidence region near the lower fault 
tip and two uplift regions on the two sides of the 
subsidence region. 

The postseismic deformations are very much affected 
by the viscosity structure in the mantle. If a low 
viscosity channel exists below the lithosphere but below 
this low viscosity channel the viscosity starts to increase, 
then the resulting time dependent deformation will be 
different from that of the low viscosity half space. In the 
low viscosity half space case, it is relatively easy for the 
material to move downwards; but if the viscosity increases 
with depth, it will then be easier for the material to move 
sideways in the low viscosity channel. Model 4B 
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demonstrates such an effect. The toodel has the same 

lithosphere and fault as in model 4A, but it has a low 

20 

viscosity channel from 80 to 160 km with viscosity 10 

22 

poise; below this depth the viscosity is assumed to be 10 
poise. The deformation at 0.0, 1.1, 3*6, and 11 years at 
the same selected locations are given in figure 4.3.5a to 
4.3.5d. The post-seismic deformation changes at 1.1, 3.6 
and 11 years at t'he same locations are given in figure 
4.3.6a to 4.3.6c. As expected, the material in the 
compressed region tends to move sideways in the low 
viscosity channel. On the free surface, the vertical 
deformation will be a smaller 'subsidence near the lower 
fault tip, because it is more difficult for the material to 
move downward; and there are two larger uplift regions 
adjacent to the subsidence region, because the material also 
flows sideways from the compressed region and material piles 
up on both sides. This phenomenon can be seen in figure 
4.3*7, where the vertical displacement change on the free 
surface are given. 

To test if the same phenomenon exists for steeply 
dipping fault, we construct two additional models. Model 4C 
and model 4D have the same fault model and lithosphere 
structure. The fault is a reverse fault dipping 60 degrees, 
the fault offset is 10 m from surface to 40 km deep and 
tapers off to zero at 45 km deep. Model 4C simulates an 
elastic lithosphere over a viscoelastic half space. The 
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elastic lithosphere is 80 km thick and has viscosity 10 

poise. Below 80 km the viscosity is 10 poise. Model 4D 

has a low viscosity channel from 80 to 180 km deep, with 
20 

viscosity 10 poise. Below 180 km deep the viscosity is 

22 

assumed to be 10 poise. The schematic diagram of the 
fault and the viscosity model are given in figure 4.3.8. 

The finite element grid used in these calculation are given 
in figure 4.3.9. 1050 degrees of freedom and 496 elements 

are used in these models. The displacement patterns for 
model 4C at time = 0.0, 5.2, and 25.7 years after the event 
are given in figures 4.3.10a to 4.3.10c, and the 
displacement changes at time = 5.2 and 25.7 years are given 
in figure 4.3.11a and 4*3. 11b, The basic pattern of the 
postseismic deformation (figure 4.3.11) is similar to that 
of the shallow dipping model 4A (figures 4.3.3a to 4.3.3c). 
Two circuits of movement are formed on the two sides of of 
the fault. This similarity indicates that the same 
mechanism is operating in both cases, namely, the material 
moves away from the compressed region below the fault and 
moves toward the fault zone from the two sides. The 
resulting vertical deformation changes for model 4C on the 
free surface are given in figure 4.3.12. As expected, there 
is subsidence near the fault and uplifts on both sides of 
the subsidence region. 

The displacement patterns for model 4D at time = 

0.0, 5.2, and 25.7 years after the event are given in figure 
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4.3.13a to 4.3.13c, and the displacement change patterns at 
time = 5.2 and 25.7 years a**e given in figures 4.3. 14a and 
4.3.14b. Here the basic mechanism for the postseismic 
movement is similar to that of model 4B. The material below 
the fault tends to flow sideways because the viscosity 
increases with depth. The vertical displacement changes at 
5.2 and 25.7 years after the event are given in figure 
4.3.15. Compared with the half space mode 4C, the 
subsidence is smaller while the uplifts are larger. 
Increasing viscosity with depth thus has evident effects on 
the deformation patterns for both shallow and steep dipping 
events . 

Thatcher and Rundle (1979) calculated the 
deformation following a thrust event in a viscoelastic 
medium. Thei»" model used an elastic layer over a 
viscoelastic half space. They found that for events that 
partially fracture the lithosphere, there is large 
subsidence near the fault zone due to viscoelastic 
relaxation after an dip slip earthquake. They also pointed 
out that the region below the fault zone is under 
compression and the material there tends to flow away, which 
will cause subsidence. They also found that buried reverse 
slip causes uplift. They proposed a model for the 
earthquake cycle in a subduction zone: The uplift after 
great thrust events is mainly due to deep aseismic slip 
while viscoelastic relaxation causes the interseismic 
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subsidence. From the examples in this section, we see that 
the large extent of subsidence caused by relaxation 
partially depends on the model. It must be relatively easy 
for the material to flow downwards to cause a large extent 
of subsidence. The deformation patterns after an earthquake 
are rather sensitive to viscosity structure. The wain cause 
of post-seismic deformation lies in the low viscosity zone: 
when the material in the low viscosity zone moves in 
response to the stress state there, the high viscosity 
regions deform accordingly in order to maintain continuity 
and equilibrium. In a subduction zone, the viscosity 
structure is much more complex than a layer over a half 
space. The complexity in viscosity structure may produce 
different postseismic deformation patterns. With realistic 
conditions for a subduction region, post-seismic uplift in a 
fault zone may be explained by viscoelastic relaxation. We 
will discuss this further in section 4.5 when we construct 
models of the 1964 Alaskan earthquake. 
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4.4 Interaction of deep aseismic slip and 
a shallow earthquake 

In this section, we study the interaction of deep 
aseismic slip and shallow earthquake along the same fault in 
a subduction region. Premonitory slow slips were observed 
before abrupt slips in rock mechanics experiments (Byerlee 
1967, Scholz et al 1972, Byerlee and Summers 1975, 

Dietereich 1978, 1979ab) . Aseismic slip on surface has been 
observed along the San Andreas fault zone (Nason 1973, King 
et al 1973, Nason and Weertman 1973). Many investigators 
believe that deep aseismic slip relieves the stress below 
the stick slip layer in San Andreas fault zone (Brace and 
Byerlee 1970, Thatcher 1975ab, Rundle and Jackson 1977, 
Turcotte et al 1979). In a subduction region, aseismic slip 
may also play a role both before and after the earthquake 
(Fitch and Scholz 1971, Scholz 1972, Kanamori 1973, Kasahara 
1975, Brown et al 1977, Fujii 1978, Fujii and Nakane 1979, 
Thatcher and Rundle 1979). It is conceivable that aseismic 
slip happens in the lower portion of the fault zone before 
an earthquake, especially if there are forces acting on the 
fault zone from the deeper part of the earth either by 
gravitational sinking or convection (Isacks and Molnar 1971, 
Toksoz et al 1973, Forsyth and Uyeda 1975, Solomon et al 
1975, Richardson 1978a). The aseismic slip loads the upper 
stick slip region, which then ruptures during an earthquake. 
Deep aseismic slip has also been used to explain the crustal 
movements after the earthquake. The possible reason for 


postseismic slip is that an earthquake loads the deeper 
portion of the fault, and the subsequent aseismic slips 
relieve the stress loaded by the earthquake. 

Figure 4.4.1 and 4.4.2 give two examples of the 
theoretical vertical deformation on the free surfac due to 
buried slips. Model 4E has a fault dipping 21.8 degrees. 

The slipping fault surface is located between the depths 40 
km and 80 km. The maximum fault slip is 1 m and the slip 
tapers off near the fault ends, as shown in figure 4.4.1. 
Model 4F has a fault dipping 45 degrees and is located 
between the depths of 40 km and 80 km, also with a maximum 
fault slip of 1 m (figure 4.4.2). The surface deformations 
in both models are mainly uplift. In the case of a normal 
fault, the vertical deformation will be mainly subsidence 
(except for the special case of vertically dipping faults). 
Thus a buried slip can be conveniently used to explain 
either uplift or subsidence. Since a buried slip cannot be 
observed directly, we have the freedom to choose its size, 
location, slip sense, and dip angle. Aseismic slip with 
properly chosen parameters can thus explain many 
observations. Discreet judgements must be used so that 
these inter pretations make physical sense. 

In the following, we construct two kinematic models 
for the interaction of deep aseismic slip with an earthquake 
along the same fault. We are interested in modelling the 


130 


* 


* 


way the aseismic slip loads the stick slip region before the 

earthquake, and the way the earthquake loads the aseismic 

slip region after the event. The finite element grid used 

in this study is given in figure 4.4.3. For the purpose of 

comparison, first let's examine the behavior of a fault in 

an elastic medium. Model 4G has elastic constants of 

12 2 

Young’s modulus 2.0 x 10 dyne/cm and Poisson's ratio 0.25 
as in other models in this chapter. The fault extends from 
surface to 80 kin depth dipping 45 degrees. The lower 40 km 
is assumed to undergo slow aseismic slip while the upper 40 
km will rupture abruptly. The time dependence of the slip 
is given in figure 4.4.4. We assume that at the time zero 
the stress in the fault zone reaches a critical level and 
the lower portion of the fault starts aseismic slip. The 
slip is assumed to be a reverse slip along the 45 degree 
fault dip. In a numerical scheme, we have to discretize the 
continuous aseismic slip, we assume the slip between the 
depths 45 km and 75 km increases from 0 to 5 m in 40 steps 
during the first 20 years. The fault slip tapers off 
between 40 and 45 km deep and between 75 and 80 km deep. At 
20 years the upper section ruptures abruptly with 5 m 
reverse slip, thus bringing equal slip displacement from the 
surface to 75 km deep, and tapers off from 75 to 80 km deep. 
The shear stress parallel to the dip direction at 6 
locations along the fault zone are given in figure 4.4.5 and 
4.4.6. Positive shear stress tends to produce reverse slip 
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on fault while negative shear stress tends to relieve the 
prestress for a thrust fault. Figure 4.4.5 shows the shear 
stress at three locations along the stick slip region. The 
aseismic slip continuously loads the stick slip region for 
20 years until the earthquake happens and the stress is 
relieved (becomes negative). Since there is no relaxation, 
the stress stays the same after 20 years. In figure 4.4.6, 
the shear stress along the aseismic slip region is given. 

The aseismic slip continuously relieves the stress in that 
region until at 20 years the earthquake happens and reloads 
the aseismic region. For this elastic medium, the load from 
the earthquake is not sufficient to bring the aseismic 
region back to zero stress (i.e. the reference stress level 
which caused the aseismic slip to happen). Actually the 
stress state afte** the earthquake is simply that of a 
uniform slip from 0 to 75 km deep and tapers off from 75 to 
80 km deep, because there is no relaxation in the model. 

Thus everywhere along the fault zone the stress is relieved. 

In model 4H we introduce viscoelastic relaxation. 

The rheological structure is given in figure 4.4.7. The 

lithosphere is assumed to be 80 km thick and has a viscosity 
25 

of 10 poise. There is a descending slab, also with the 
25 

viscosity 10 poise. Outside the descending slab, the 

20 

viscosity is 10 poise from 80 to 160 km depth, and below 
1 6 Q km the viscosity is assumed to be 10 poise. We also 
assume there is a low viscosity wedge above the shear zone 
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between the descending and overiding lithospheres, and the 

20 

viscosity is assumed to be 10 poise there. The resulting 
shear stress parallel to the dip direction at the same 
locations are given in figure 4.4.8 and 4.4.9. In the stick 
slip region, the results are similar to those of the elastic 
model (figure 4.4.8), except that the magnitude of loading 
from the aseismic slip is slightly reduced, and after the 
earthquake the stress slowly recovers towards zero. In the 
aseismic slip zone, where we assumed low viscosity, the 
shear stress is quite different from the elastic case. In 
the first 20 years, the effect of viscoelastic interaction 
is to recover the prestress. The end result is that the 
stress relief at 20 years is considerably less than that of 
the elastic case. The earthquake can then load the upper 
part of the aseismic slip zone back to positive stress (i.e. 
stress larger than the reference stress which started the 
aseismic slip) . The behavior after the earthquake is that 
shear stress relaxes towards zero. 

Admittedly this simple model is far from perfect. 

For Example, it is probably not reasonable to assume that 
the aseismic slip starts simultaneously along the whole 
aseismic region. Maybe it should start from one end and 
propagate to the other end. We will not attempt to model 
this mode (and infinitely many other modes) of slips in this 
thesis. However, this simple numerical experiment does tell 
us something about the interaction between the earthquake 
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and the aseismic slip: the aseismic slip stresses the stick 
slip zone before the earthquake, the earthquake in turn 
loads the aseismic slip zone. The stress caused by the 
earthquake on the aseismic zone usually decreases with 
distance from the source. Viscoelastic interaction with the 
aseismic slip has a stress recovery effect in the aseismic 
slip zone. Combination of the loads from the earthquake and 
the stress recovery due to viscoelastic relaxation can 
sometimes bring the region near the earthquake fault tip to 
the stress level which started the aseismic slip. 
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4.5 Modelling of the 1964 Alaskan earthquake 

4.5.1 General descr itptions 

The coseismic deformations of the 1964 Alaska 

earthquake are very complex. Many subsidiary faults were 

activated during the earthquake and contributed locally to 

the deformations (Plafker 1969). On the other hand, 

repeated leveling data showed a rather smooth profile over a 

long levelling route (figure 4.2.5); this suggested the 

nonlocal nature of the time dependent deformation. This 

deformation was probably associated with the main fault 

rather than the subsidiary local effects. We set up a model 

to calculate the relaxation associated with the main fault 

of the 1964 earthquake. The observed data indicated that 

the variation of uplift in the direction perpend icular to 

the leveling route was small (figure 4.2.5). Thus we can 

use two dimensional plane strain models. All the models in 

this section have the same finite element grid, fault 

parameters, and instataneous response. The finite element 

grid has 1412 degrees of freedom and 651 elements (see 

figure 4.5.1). The model region is 2200 km long and 1000 km 

deep. As usual we assume that the material in the models is 

elastic in bulk and Maxwellian viscoelastic in distortion. 

The elastic parameters are assumed to be Young’s modulus 2.0 
12 2 

x 10 ' dyne/cm and Poisson's "'atio 0.25. The differences 

between the models will be due to different viscosities. 
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The purpose of this study is not to fit the geodetic 
data in detail with a complicated model. We are interested 
in using the features in the data to bring out the physical 
processes occurring after the earthquake. The models are 
constructed according to our knowledge of subduction 
regions. Starting with a simple layered model, one by one 
we add the features observed in a subduction region to the 
simple model. We then compare the model results with the 
uplift data. Hopefully with this process we can gain some 
insights into the physics of the time dependent movements in 
a subduction region. 

4.5.2 Fault model 

The schematic model of the fault is shown in figure 
4.5.1. The values for the maximum fault slip (12 m) and the 
fault width (203 km) are similar to those of the optimal 
model of Hastie and Savage (1970). However, the dip angle 
of the fault plane in Hastie and Savage’s model is only 3-7 
degrees, while we used a dip angle of 10 degrees. Our value 
is closer to the values determined by seismic studies (20 
degrees from surface wave analysis by Kanamori, 5 to 15 
degrees from body wave analysis by Stauder and Bollinger) . 

We also assume that the fault slip tapers off linearly to 
zero near where the fault ends instead of ending abruptly 
with a uniform slip. As in the optimal model by Hastie and 
Savage, we assume that there is no surface rupture. The 


topmost point of the fault is 5 km deep. There are few data 
near the trench, and it is not certain if the main fault 
reached the surface. However, the detailed configuration of 
the fault near the upper fault tip should have little effect 
on the deformation from Whittier to Anchorage, which are 
located near the lower tip of the fault some 200 km away. 

We chose a buried fault partly for computational economy. 

In order to maintain the proper aspect ratio of the elements 
in the wedge area above a shallow dipping fault, many 
elements are needed near the intersection of the fault and 
the free surface if we use quadrilateral elements. The grid 
structure is simpler for a buried fault. In the coordinate 
system of the finite element grid, the fault is located 
between the horizontal coordinates -200 km and 0 km, and 
between 5 and MO km of depth. 

4.5.3 Coseismic response 

All the models have the same elastic parameters, 
thus they have the same coseismic responses. Figure 4.5.2 
shows the coseismic vertical deformation on the free surface 
from the model calculation. The general behavior is 
subsidence in the inland side and uplift toward the trench 
on the overiding lithosphere, and little deformation on the 
footwall side. This is expected for a shallow dipping 
dislocation (Plafker 1972, Mansinha and Smylie 1971). The 
vertical deformation does not have the spectacular uplift 
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associated with steep dipping subsidiary faults as in the 
Hastie and Savage model, but the general features are 
similar. The maximum uplift is located about 30 km from the 
fault tip. The rather rapid increase of uplift there is due 
to the material pile up near the buried fault tip. The same 
fault model has also been used in section 4.3 to study the 
deformation patterns of the shallow dipping event. The 
displacement vectors at selected nodes were given in figure 
4.3.2a of section 4.3. 

The town of Whittier is located near the maximum 
coseismic subsidence along the leveling. In the model 
result, this point is near the horizontal coordinate -180 
km. We will use the segment from -280 km to -180 km as the 
leveling route when comparing the model results with the 
uplift data. 

4.5.4 The layered models 

For comparison purposes, we calculated the time 

dependent response of two simple layered models. Model 41 

is identical to model 4B in section 4.3. We assumed that 

25 

the lithosphere is 80 km thick and has a viscosity of 10 

poise. Below the lithosphere is a layer of low viscosity 

10 extending from 80 to 160 km deep. Below 160 km the 
22 

viscosity is 10 poise. The vertical displacement changes 
at 1.1, 3-6 and 11 years after the event are given in 
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figure 4.5.3* Superimposed on this figure is the uplift 
data of surveys in 1965, 1968 and 1975. Theoretical results 
are subsidence near the lc-wer tip of the fault and uplift on 
both sides of the subsidence area. It is interesting that 
the uplift on the inland side has a magnitude of about 60 cm 
in 11 years, comparable to the observed uplift data. On the 
whole though, the match between the model result and the 
observation is not satisfactory. The position of the inland 
uplift in the model is located more to the inland side. 
Nevertheless, this simple model illustrates an important 
point: for reasonable values of viscosity, great earthquakes 
like the 1964 Alaska earthquake can produce measurable 
deformation years after the event. In the case of the 
Alaska earthquake, the amount of uplift in the model is 
comparable to the leveling result. This suggests that it is 
important to include the effects of viscoelastic relaxation 
in the analysis of time dependent geodetic measurements. 

Model 4J is identical to model 4A in section 4.3, it 

simulates an elastic lithosphere over a viscoelastic half 

space. This model differs from model 41 only in that the 

viscosity below 160 km is assumed to be 10^ poise instead 
22 

of 10 poise. As discussed in section 4.3, this would 
cause a larger subsidence zone because it is easier for the 
material to move downwards. The result of the vertical 
displacement change is given in figure 4.5.4. The amplitude 
and the extent of the subsidence are larger than those in 
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the layered ease model 41; end the uplift in the inland side 
is reduced (notice the scale change in the figure). The 
results of this model are still not compatible with the 
observations, however, the uplifts are still of measurable 
magnitude . 

4.5.5 Model with descending lithosphere 

In model 4K, we introduce a short descending 
lithosphere (see figure 4,5.5). Alaska is a part of the 
circum-Pacif ic seismic belt. Seismicity studies indicated 
that intermediate earthquakes exist in this region 
(Gutenberg and Richter 1954, Isacks and Barazangi 1977, Van 
Wormer et al 1974); therefore the existence of a descending 
slab is expected from plate tectonics theory. We assumed 
that a short descending slab exists in model 4K. The 
descending lithosphere is assumed to be 80 km thick, and the 
maximum depth of the descending slab is 1 60 km. We assume 
that the viscosity within the descending slab is the same as 
the lithosphere above 80 km. The viscosity in the other 
area is the same as in the layered model (model 41). The 
result for vertical deformation is given in figure 4.5.6. 
Compared with the result of model 41 (figure 4.5.3), both 
the uplifts in the inland and the subsidence are reduced. 

The physical reason for this difference can be explained in 
the following way: Since the descending lithosphere blocks 

the low viscosity channel, it is more difficult for the 
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material under the fault zone to move to the inland side. 
Thus the uplift in the inland side is reduced. The 
subsidence near the fault is also reduced, because less time 
dependent movement can be induced due to the high viscosity 
descending lithosphere. As seen in figure 4.5.6. This 
model result also does not explain the uplift data. 

4.5.6 Model with descending slab and low viscosity wedge 

The wedge areas above the descending slabs often are 

regions of high heat flow (Watanabe et al 1977) and high 

seismic wave attenuation (Oliver and Isacks 1967, Molnar and 

Oliver 1971, Utsu 1971, Barazangi et al 1975). Shear 

heating between the descending and overiding lithospheres 

may have raised the temperature in this region (Toksoz et al 

1973, Hsui and Toksoz 1979). Since viscosity depends 

strongly on temperature, it is conceivable that the region 

above the shear zone can have a rather low viscosity. In 

addition, if strain rate depends on higher power of stress, 

the effective viscosity near the fault zone will be lower 

because the stress is higher there. In model 4L'we- assume 

that the low viscosity extends to a thin layer above the 

boundary between the descending and overiding lithospheric 

plates, as shown in figure 4.5.7. We assumed that the 

1 9 

viscosity in this low viscosity zone is 2.0x10 poise. 

Other aspects of this model are the same as those for model 
4K . 
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The deformation changes on the free surface at 1.1, 
3.6 and 11 years after the event are given in figure 4.5.8, 
Superimposed on this figure a**e the leveling data from 
1965,1968 and 1975 surveys. The features of the model 
results and the observations match very well. The basic 
features of the theoretical displacement profile are still 
two uplift areas separated by a subsidence region near the 
fault tip, as in the case of layered model 41 (cf. figure 
4.5.3), But the uplift region on the inland side is 
narrower than that in the layered model. The position of 
the uplift has also shifted towards the trench side. The 
points of maximum uplift migrate towards the trench with 
time. This feature also exists in the uplift data, as 
pointed out by Brown et al (1977), Figure 4.5.9 shows the 
enlarged figure of the model results superimposed on the 
observed uplift data . The match at 11 years after the 
event is excellent. Such a good match is rare in the 
modelling of geodetic data, and probably is fortuitous, 
because the true structure of the earth is certainly more 
complex than the simple model. Nevertheless, it is 
gratifying to find that a model constructed from our 
knowledge about the subduction region can produce the 
features in the data, and gives better fit than other less 
realistic models. The model uplift at 1.1 and 3.6 years 
after the event is smaller than the 1965 and 1 96 8 survey 
data. This is not surprising though: We have assumed that 

the earthquake happens instataneously . For the first 
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several years, the crustal movements associated with the 
aftershocks and deep slips probably contributed to the 
observed deformation. The relaxation effect became dominant 
only after the phase of rapid adjustments like aftershocks 
had decayed. Notice that the discrepancies between data and 
theoretical results is larger near the fault end (which is 
located at -200 km in the horizontal coordinate) than near 
the inland side for the early years. Some extension of the 
earthquake fault may have occurred through aseismic slip. 

The buried aseismic slip in general will cause uplift above 
it, as we have seen in section 4.4. It is also reasonable 
that the aseismic slip is larger near the fault end, where 
the stress caused by the earthquake is largest. 


The displacement vectors at selected locations at 
tia;e = 0.0, 1.1, 3.6, and 11 years after the event are given 
in figure 4.5.10a to 4.5. lOd. The general pattern is still 
that of a double couple source. However, the downward 
displacements under the lower fault tip in general decrease 
rather than increase with time, indicating that the material 
there moves upwards rather than downwards after the event. 
The direction of the material movements in the compressed 
region is thus opposite to those in the cases of the half 
space and the layered model (model 4A and 4B in section 
4.3). This behavior is more clearly seen in figure 4.5.11a 
to 4.5.11c, where the displacement changes at 1.1, 3.6, and 
11 years are given. This deformation pattern is probably 
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due to the interaction of the low viscosity zone and the 
high viscosity descending slab. Because the high viscosity 
descending lithosphere underneath the low viscosity region 
serves as a guide for the material movements, the material 
in the low viscosity zone above the descending slab cannot 
move downwards into the descending slab; it will move either 
downdip or updip. To maintain continuity and equilibrium 
the materials above and below this low viscosity region will 
move towards it. This created the complex deformation 
pattern in figure 4.5.11. To show the displacement change 
pattern more clearly, we reduced the scale of displacement 
and replotted figure 4.5.11c in figure 4.5. lid. Notice that 
in the downdip extension region of the fault, the 
deformation pattern is very similar to that of a reverse 
fault along the segment indicated by the dashed line. This 
is hardly surprising, as this configuration tends to relieve 
the earthquake load in that region. 

4.5.6 The effects of slab geometry 

In model 4K and 4L, we have kept the shape of the 
descending slabs simple. The real situation may be more 
complicated. Seismicity profiles in this region are rather 
diffused (Van Wormer et al 1974). Van Wormer et al (1974) 
suggested that there are two blocks in South central Alaska. 
The boundary between the two blocks extends from the Yetna 
River through Prince William Sound, perhaps as far as the 
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continental shelf break near Yakutat Bay. The leveling 
route discussed in this section and the epicenter of the 
1964 Alaskan earthquake are located near this boundary. The 
subduction direction north of this block boundary is about 
N40W, while to the south of the boundary it is about N77W. 

We cannot take this complexity into account in our two 
dimensional models. However, this difference was inferred 
from the seismicity in the intermediate depth. If any split 
of the descending slab does exist, it should exist only in 
the deeper part of the slab. In the shallow part of the 
descending slab it should be continuous. In model 4M we 
investigate the effects of the deep structure of the 
descending slab on the relaxation. The model is given in 
figure 4.5.12. Model 4M is the same as model 4L except for 
the addition of a deeper descending slab. The result of the 
surface deformation is given in figure 4.5.13. The uplift 
data are also shown on the figure. The effect of the deeper 
part of the descending lithosphere on the surface 
deformation is very small. On the other hand, this also 
means that it is not very useful to use surface deformation 
to probe the deep structure of a descending lithosphere. 


4.5.7 The effect of viscosity 


The viscosity in the earth is still a controversial 
issue. Since viscosity depends strongly on temperature and 
temperature anomalies exist in subduction regions, the 
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viscosity in a subduction region may be quite different from 

that in other areas. The value of viscosity in a subduction 

region is largely speculative. Several investigators have 

pointed out that great thrust earthquakes in a subduction 

zone may provide an opportunity to investigate the viscosity 

there (Nur and Mavko 1974, Smith 1974, Kasahara 1975, 

Thatcher and Rundle 1979). We have found for the case of 

the Alaskan earthquake, the leveling data put some 

constraint on the possible models. The models with a 

descending slab and low viscosity wedge above the descending 

slab can reproduce the features in the data, while the 

simple layer models and the model without a low viscosity 

wedge produced less satifactory results. To see the 

sensitivity of surface deformation on the viscosity, we 

constructed one additional model with different viscosity in 

the low viscosity channel. Model 4N is similar to model 4L, 

1 9 

except the low viscosity below 80 km deep is 5 x 10 poise 
20 

instead of 10 poise. The result of the vertical 
displacement changes are given in figure 4.5.14, and the 
uplift data are superimposed on the model result. This 
model gives a less satisfactory fit to the uplift data than 
model 4L. The result has a broader extent of uplift 
compared with the data, and the uplift amplitude in the 
model is larger. This illustrates that the surface 
deformation is sensitive to the viscosity structure, and 
geodetic measurements can be used to constrain the possible 
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4.6 Discussion 

In this chapter, we will look at the data and models of 
the time dependent deformations in subduction zones. In 
general, a region with a long subduction history like 
Hokkaido and North Honshu has a rather complex structure. 

The origin of the tectonic features are not fully understood 
at the present time. The geodetic data also cannot be 
easily explained by a simple underthrust and rebound cycle. 
On the other hand, in Southwest Japan and Alaska, geodetic 
data seem to suggest a simple picture of an earthquake 
cycle: Th ; e coastal region is dragged down by the descending 

lithosphere before the great shallow earthquake, and it 
rebounds upward following the earthquake. We looked into 
some basic physical processes in this cycle, In section 
4.3, we noticed that the time dependent deformations are 
very sensitive to the viscosity distribution. The movements 
induced by the earthquake originate primarily in the low 
viscosity zone. The high viscosity regions deform in 
response to the movements in the low viscosity region to 
maintain continuity and equilibrium. This sensitivity of 
deformation to viscosity makes geodetic data a valuable tool 
for probing the viscosity structure of a fault zone. In 
section 4.4, we looked into the interaction of preseismic 
slow slip and the earthquake. We found that deep preseismic 
slip stresses the earthquake fault zone and the earthquake 
in turn will load the deep aseismic slip zone. Viscoelastic 
interaction with the aseismic slip causes a stress recovery 
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effect in the aseismic slip region. Due to the load of the 
earthquake and the stress recovery effect, it is possible 
that afterslip occurs near the earthquake fault tip. In 
section 4.5, we constructed a series of models for the 1964 
Alaskan earthquake and compared the model results with the 
time dependent geodetic data. The most important finding in 
this modelling is that with a realistic viscosity structure 
of a subduction zone, viscoelastic relaxation can explain 
the uplift data. In the modelling of a subduction zone, it 
is essential to include a high viscosity descending slab and 
a low viscosity region above the shear zone between the 
descending slab and the overiding lithosphere. Simpler 
model results can not fit the uplift data. Since surface 
deformations are sensitive to the viscosity in the fault 
zone, geodetic data can put stringent constraints on the 
possible viscosity models. As the viscosity is vitally 
relevent to stress accumulation and earthquake occurrence, 
the modelling should be used with the high quality data from 
space techniques to investigate earth structure. 

As we have mentioned, the purpose of this study is 
not to fit the Alaskan data in detail with some complicated 
model. Rather, we are interested in bringing out the 
important physical processes for an observed time dependent 
phenomenon, and investigating the approach we can follow in 
the future when more data become available. In the process 
of finding earth structure f^om observed geodetic data, 
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there is the problem of uniqueness of the solution. With 
only one set of leveling data, there are other possible 
explanations. Brown el al (1977) examined the possible 
processes that could explain the uplift following the 1964 
Alaskan earthquake. The possible processes were (1) creep 
on a buried fault, (2) viscoelastic rebound, (3) strain 
accumulation due to lithosphere loading, (4) magma 
intrusion, and (5) dilatancy effects. Brown et al (1977) 
concluded that magma intrusion and dilatancy effects were 
not the probable explanations for the observed phenomenon: 
there was no indication of magma activities in the region; 
and dilatancy predicts subsidence rather than uplift, 
contrary to what was observed. They thought that strain 
accumulation due to lithospheric loading and viscoelastic 
rebound might play a part in the uplifting process, but they 
did not think these two processes were important. Strain 
accumulation is too slow a process to cause such a rapid 
uplift. On the viscoelastic rebound, they quoted the 
results of Nur and Mavko (1974) to estimate the effects of 
viscoelastic relaxation. The magnitude of theoretical 
uplift was found different from the observed data. This is 
not surprising since a layered model does not reflect the 
structural features of a subduction zone. In addition, the 
solution of Nur and Mavko is a doubtful one; it has been 
criticized by Smith (1974) and Thatcher and Rundle (1979). 
Brown et al (1977) also doubted if the viscosity in the 
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fault zone could be as low as 5 x 10 poise, as implied in 
the calculations by Nur and Mavko. There is no definite 
evidence to argue against or for this viscosity value in a 
fault zone. We are not aware of any heat flow measurements 
near the leveling route. There are volcanos about 100 km 
away from Anchorage. On the other hand, Barazangi et al 
(1975) studied the body wave attenuation in the Aleutian arc 
about 1000 km away from this region and found no evidence of 
a high attenuation zone there. We think the problem of the 
viscosity value in a fault zone is still open. However, 
when modelling a subduction zone, we should not be bound to 
use the viscosity value derived from the post glacial 
rebound data in continental areas. The shear heating may 
reduce the viscosity near the fault zone. In addition, the 
exact form of a constitutive relation between stress and 
strain is not known. If strain rate is proportional to a 
certain power of stress, the effective viscosity in a fault 
zone will be lowered due to the high stress there. 

The favored explanation of the uplift data by Brown et 
al (1977) is the deep aseismic slip on an extension of the 
fault plane which ruptured during the 1964 earthquake. A 
2.3 m slip dipping 25 degrees on a buried fault segment 
about 50 km wide near Anchorage gives a satisfactory fit to 
the uplift at 1975 (Brown et al 1977). In order to expalin 
the migration of maximum uplifts f^om the inland toward the 
trench, Brown et al suggested the following process: In 
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south central Alaska, the megathrust zone between the North 
American plate and the underthrusting Pacific plate can be 
divided into two major segments. The 1964 Alaskan earthquake 
ruptured the upper section of the lithosphere where the 
material is brittle. The lower segment of lithosphere is 
dominated by viscous creep. The segment of the fault which 
is inferred to be responsible for the postseismic movement 
is in a transition zone between the pure stick slip and the 
pure creep behavior. The rupture of the brittle upper 
segment loads the transition zone, which then proceeds to 
fail via viscous creep. If the time constant of the creep 
mechanism decreases with depth along the fault, the effects 
of the slip on the lower portions will be felt first and 
peak out before those of the upper portions. Smaller time 
constants imply smaller effective viscosities, which may be 
associated with the geothermal gradient in the area. The 
lower portions, being hotter, should have smaller time 
constant of creep. 

In the process proposed by Brown et al , it is 
assumed that a zone of low viscosity exists in the downdip 
extension region of the earthquake fault. The relaxation 
time of the region should be less than ten years so that the 
effects of migration of maximum uplifts can be seen within 
one decade after the event. The existence of such a low 
viscosity region is essential in the interpr etations by 
Brown et al : If such a low viscosity zone does not exist, 
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the aseismic slip should occur near the earthquake fault tip 
rather than further downdip, since the load caused by the 
earthquake decreases rapidly with distance and the 
frictional strength of the fault increases with depth. It 
would then be difficult to explain the gap between the 
aseismic slip zone and the earthquake fault; it would also 
be difficult to explain why the slip went updip. If such a 
low viscosity zone exists, the viscoelastic relaxation 
effect after the earthquake should be evaluated. In section 
4.5 we see that the deformation pattern due to relaxation in 
the downdip extension region of the fault is that of a 
reversed slip (cf. figure 4.5. lid). In a sense, the 
relaxation has produced some kind of aseismic slip, except 
that the slip does not occur d iscontinuously across a 
surface, and the deformation pattern changes continuously 
but rapidly in a zone. However, in the relaxation model, 
the trenchward migration of maximum uplift comes about as a 
natural result of the structure of the subduction zone. 

Brown et al assumed the existence of a low viscosity zone 
but did not calculate the effect of the viscoelastic 
relaxation, instead they imposed aseismic slips to fit the 
data. In a sense, the relaxation models give a more 
complete description of the physical process. We don't have 
to imposed the aseismic slip in an ad hoc fashion, and the 
viscosity model used in the calculation is a realistic model 
of the subduction zone. The roles of a descending slab and 
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a low viscosity zone above the descending slab are essential 
in the relaxation model; The descending slab guides the flow 
direction in the low viscosity zone. Had there been no 
desending slab, the vertical deformation on the surface 
would be dominated by subsidence near the lower fault tip, 
because the material movement will be mainly downwards 
there, as in the layered models in section 4.5 

Thus in some sense the relaxation model contains the 
aseismic slip model. However, the relaxation model is 
constructed with our knowledge of the subduction zone; the 
comparison of geodetic data and model-s thus provides us with 
the information on the viscosity structure in a subduction 
zone. If our interpretation is correct, the vertical 
deformation from Whittier to Anchorage should still be 
uplift in the next decade but with decreased uplift rate. 
Figure 4.6.1 gives the prediction by using the results of 
model 4L. The theoretical uplifts at 1.1, 3-6, 11, 13.2, 
19.1 and 27.5 years after the event are shown in figure 
4.6.1. It must be clarified that even if this prediction 
comes true, by using only one leveling route, it is not 
possible to acertain the viscosity structure of Alaska. We 
have few constraints on the trench side. However, the 
general features of the uplift data can be explained by 
viscoelastic relaxation in a fault zone with a descending 
slab, a low viscosity layer below the lithosphere, and a low 
viscosity zone above the shear zone. These are probably the 
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basic elements of a subduction zone structure. For more 
detailed structures of subduction regions, more geodetic 
measurements are needed* 

Addendum : 

A paper (Wahr and Wyss 1980) appeared after this thesis 

had been drafted, Wahr and Wyss (1980) modelled the post- 

seismic uplift of the 1964 Alaskan earthquake with a two 

dimensional finite element model. Their model consists of 

an elastic half space and a viscoelastic inclusion. The 

viscoelastic inclusion models the low viscosity wedge above 

the shear zone. They found that the uplift from Whittier to 

Anchorage at 1975 can be fit by a model having an inclusion 

19 19 

viscosity 1.2 x 10 to 2.2 x 10 poise. This viscosity 
value is similar to our results in this chapter. However, 
the maxima of uplift migrate away from the trench in their 
model, in the opposite direction as the data indicated. 
Models 4L and 4M in this chapter included more features of 
the subduction region and the migration direction of maximum 
uplift conformed to what was observed. 
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Figure Captions - Chapter 4 

Figure 4.2.1. Tonankai and Nankaido earthquakes. The 

aftershock area 1 day after the main shock is shown by 
solid curves and that after one month by dashed curves. 
The solid polygon shows the projection of the 
theoretical fault plane from Fitch and Scholz (1971). 
The leveling result from points a to b is given in 
figure 4.2.2. Point b is Muroto Promontory, which is 
the closest point to the trench along the leveling 
route. (Reproduced from Kanamori 1973). 

Figure 4.2.2 Pre-, co- and post-seismic vertical 

deformation along the route from a to b in figure 
4.2.1. (Reproduced from Kanamori 1973). 

Figure 4.2.3 Observed and calculated vertical 

displacements for the 1964 Alaska earthquake along a 
cross section through Montague and Middleton Islands. 
(Reproduced from Hastie and Savage 1970). 

Figure 4.2.4 Simplified tectonic map of South Central 

Alaska. The epicenter of the 1964 Alaska earthquake is 
denoted by a circled star. The route of the leveling 
profile is shown as a light dashed line. (Reproduce 
f^om B^own et al 1971). 

Figure 4.2.5 Profile of the elevation change along the 

leveling line shown in figure 4.2.4. (Reproduced from 
Brown et al 1977). 
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Figure 4.3.1 (a) The finite element grid used for model 4A 
and 4B. The grid region is 2,200 km long and 1,000 km 
deep. It has 1412 degrees of freedom and 651 elements, 
(b) The fault model used for model 4A and 4B. The 
fault is located between -200 km and 0 km in the 
horizontal coordinate and between 5 and 40 km in depth 
in the finite element grid coordinate system. The 
fault has 12 m reverse slip along the solid line 
segment and the slip tapers off near the ends along the 
segments of dashed lines. 

Figure 4.3*2 (a) The displacement vectors at selected 
locations immediately after the event for model 4A. 

The small circles indicate the locations where the 
calculations are made. The small line segments 
stemming from the circles indicate the displacements at 
the locations. A scale of the displacement is given in 
the lower left corner and a scale of the map is given 
in the lower right corner. The thick line indicates 
the location of the earthquake fault. (b) The same 
plot at 1.1 years after the event. (c) The same 
plot at 3.6 years after the event. (d) The same plot 
at 11.0 years after the event. 

Figure 4.3*3 (a) The changes of displacement vectors at 
selected locations at 1.1 years after the event for 
model 4A . The small circles indicate the locations 
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whe^e the calculations are made. The small line 
segments stemming from the circles indicate the changes 
of displacements at the locations. A scale of the 
displacement is given in the lower left corner and a 
scale of the map is given in the lower right corner. 

The thick line indicates the location of the earthquake 
fault. (b) The same plot at 3.6 years after the 
event. (c) The same plot at 11.0 years after the 
event . 


Figure 4.3.4 The vertical displacement 
surface at 1.1, 3.6 and 11.0 years 
model 4A 


changes on the free 
after the event for 


Figure 4.3.5 (a) The displacement vectors at selected 
locations immediately after the event for model 4B. 

The small circles indicate the locations where the 
calculations are made. The small line segments 
stemming from the circles indicate the displacements at 
the locations. A scale of the displacement is given in 
the lower left corner and a scale of the map is given 
in the lower right corner. The thick line indicates 
the location of the earthquake fault. (b) The same 


t 

at 

1.1 years a 

fter 

the event. 

(c) 

The same 

t 

at 

3.6 years a 

f ter 

the event. 

(d) 

The same plot 

1 1 

. 0 

years afte-' 

the 

event. 




Figure 4.3*6 (a) The changes of displacement vectors at 
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selected locations at 1.1 years after the event for 
model 4B. The small circles indicate the locations 
where the calculations are made. The small line 
segments stemming from the circles indicate the changes 
of displacement at the locations. A scale of the 
displacement is given in the lower left corner and a 
scale of the map is given in the lower right corner. 

The thick line indicates the location of the earthquake 
fault. (b) The same plot at 3.6 years after the 
event. (c) The same plot at 11.0 years after the 
event . 


Figure 4.3.7 The vertical displacement 
surface at 1.1, 3.6 and 11.0 years 
model 4B Figure 


changes on free 
after the event for 


Figure 4.3.8 The schematic diagram of the fault and the 
viscosity structure for model 4C . The fault has 
uniform slip of 10 m from 0 to 40 km deep and tapers 
off to zero at 45 km deep. The numbers with exponents 
are viscosity values in poise. 

Figure 4.3.9 Finite element grid for models 4C and 4D. 

Figure 4.3*10 (a) The displacement vectors at selected 
locations immediately after the event for model 4C . 

The small circles indicate the locations where the 
calculations are made. The small line segments 
stemming from the circles indicate the displacements at 


the locations. A scale of the displacement is given in 
the lower left corner and a scale of the map is given 
in the lower right corner. The thick line indicates 
the location of the earthquake fault. (b) The same 
plot at 5.2 years after the event. (c) The same 
plot at 25.7 years after the event. 

Figure 4.3.11 (a) The changes of displacement vectors at 
selected locations at 5.2 years after the event for 
model MC. The small circles indicate the locations 
where the calculations are made. The small line 
segments stemming from' the circles indicate the changes 
of displacements at the locations. A scale of the 
displacement is given in the lower left corner and a 
scale of the map is given in the lower right corner. 

The thick line indicates the location of the earthquake 
fault. (b) The same plot at 25.7 years after the 
event . 

Figure 4.3.12 Vertical displacement change on the free 

surface at 0.0, 5.2 and 25.7 years after the event for 
model 4C. 

Figure 4.3.13 (a) The displacement vectors at selected 
locations immediately after the event for model 4D. 

The small circles indicate the locations where the 
calculations are made. The small line segments 
stemming f^om the circles indicate the displacements at 
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the locations. A scale of the displacement is given in 
the lower left corner and a scale of the map is given 
in the lower right corner. The thick line indicates 
the location of the earthquake fault. (b) The same 
plot at 5.2 years after the event, (c) The same plot at 
25.7 years after the event. 

Figure 4.3.14 (a) The changes of displacement vectors at 
selected locations at 5.2 years after the event for 
model 4D. The small circles indicate the locations 
where the calculations are made. The small line 
segments stemming from the circles indicate the changes 
of displacements at the locations. A scale of the 
displacement is given in the lower left corner and a 
scale of the map is given in the lower right corner. 

The thick line indicates the location of the earthquake 
fault. (b) The same plot at 25.7 years after the 
ev ent . 

Figure 4.3.15 Vertical displacement change on the free 

surface at 0.0, 5.2 and 25.7 years after the event for 
model 4D. 

Figure 4.4.1 The vertical displacement on the f^ee surface 
for a buried fault model (model 4E). The fault model 
- is given in the lower figure. The fault is located 
, ; between 40 and 80 km deep. The dip angle is 21.8 

degrees. The maximum slip is 1 m and the fault slip 
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tapers off along the dotted line near the ends of the 
fault , 

Figure 4.4.2 The vertical displacement on the free surface 
for a buried fault model (model 4F) . The fault model 
is given in the lower figure. The fault is located 
between 40 and 80 km deep. The dip angle is 45 
degrees. The maximum slip is 1 m and the fault slip 
tapers off along the dotted section near the ends of 
the fault. 

Figure 4.4.3 The finite element grid for model 4G and 4H. 

Figure 4.4.4 (a) The time dependence of the slip between 
45 and 75 km deep of the fault in model 4G and 4H. The 
slip tapers off to zero between 45 and 40 km deep and 
between 75 and 80 km deep, (b) The time dependence of 
slip along the upper 40 km of the fault in model 4G and 
4H . 

Figure 4.4.5 The shear stress parallel to 45 degree dip at 
three locations in the stick slip region for model 4G . 
The locations are indicated in the inset. Positive 
stress implies that stress is accumulating and negative 
stress implies that stress is relieved from the 
reference stress state. 


Figure 4.4.6 The shear stress parallel to 45 degree dip at 


162 


three locations in the aseismic slip region for model 
4G . The locations are indicated in the inset. 

Figure 4.4.7 Schematic diagram for model 4H. The hatched 

20 

area has a low viscosity 10 poise. The lithosphere 

25 

and the descending slab has a viscosity of 10 poise. 
The fault dip is 45 degrees and extends f^om 0 to 80 
km. The lower 40 km slips aseismically ; the upper 40 
km ruptures abruptly during the earthquake. 

Figure 4.4.8 The shear stress parallel to 45 degree dip at 
three locations in the stick slip region for model 4H . 
The locations are indicated in the inset. 

Figure 4.4.9 The shear stress parallel to 45 degree dip at 
three locations in the aseismic slip region for model 
4H. The locations are indicated in the inset. 

Figure 4.5.1 (a) The finite element grid used for studying 
the 1964 Alaska earthquake. The grid region is 2,200 
km long and 1,000 km deep. It has 1412 degrees of 
freedom and 651 elements. (b) The fault model used 
for studying the 1964 Alaska earthquake. The fault is 
located between -200 km and 0 km in horizontal 
coordinate and between 5 and 40 km in depth in the 
finite element grid coordinate system. The fault has 
12 m reverse slip along the solid line segment and the 
slip tapers off near the ends along the segments of the 
dashed lines. 
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Figure 4.5.2 The coseismic vertical displacement on the 
free surface for the fault model in figure 4.5.1 

Figure 4.5.3 The vertical displacement change on the free 
surface for model 41. The observed uplift data f'-om 
Brown et al are also superimposed on the figure. The 
model results are at 1.1, 3.6 and 11 years after the 
event for curves 1, 2 and 3. The observed uplifts at 
1965, 1968 and 1975 are denoted by crosses, circles and 
stars . 

Figure 4.5.4 The vertical displacement change on the free 
surface for layered model 4J. The observed uplift data 
from Brown et al are also superimposed on the figure. 
The model results are at 1.1, 3.6 and 11 years after 
the event. The observed uplifts are from 1965, 1968 
and 1975. Notice the vertical scale change from figure 
4.5.3. 

Figure 4.5.5 The viscosity distribution for model 4K. 

There is a short descending lithosphere to 160 km deep. 

25 ' 

The descending lithosphere has a viscosity 10 poise. 
Other parameters are the same as model 41. 

Figure 4.5.6 The vertical displacement change on the free 
surface for model 4K. The observed uplift data from 
B»"Own et al are also superimposed on the figure. The 
model results a^e at 1.1, 3.6 and 11 years after the 
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event. The observed uplifts are at 1965, 1968 and 
1975. 

Figure 4.5.7 The viscosity distribution for model 4L. The 

descending slab is the same as in model 4K, but there 

is a low viscosity zone above the shear zone between 

the descending and overiding lithospheres. The 

1 Q 

viscosity is assumed to be 2.x 10 7 poise there. 

Figure 4.5.8 The vertical displacement change on free 

surface for model 4L. The observed uplift data from 
Brown et al are also superimposed on the figure. The 
model results are at 1.1, 3.6 and 11 years after the 
event. The observed uplifts are at 1965, 1968 and 
1975. 

Figure 4.5.9 The vertical displacement change from -280 km 
to -180 km for model 4L. Results are enlarged from 
figure 4.5.8. The uplift f^om Whittier to Anchorage is 
superimposed on the figure. 

Figure 4.5.10 (a) The displacement vectors at selected 
locations immediately after the event for model 4L. 

The small circles indicate the locations where the 
calculations are made. The small line segments 
stemming from the circles indicate the displacements at 
the locations. A scale of the displacement is given in 
the lower left corner and a scale of the map is given 
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in the lower right corner. The thick line indicates 
the location of the earthquake fault. The boundaries of 
low viscosity zones and the descending slab are also 
indicated on the plot, (b) The same plot at 1.1 years 
after the event. (c) The same plot at 3.6 years 
after the event, (d) The same plot at 11.0 years 
after the event. 

Figure 4.5.11 (a) The changes of displacement vectors at 
selected locations at 1.1 years after the event for 
model 4L. The small circles indicate the locations 
where the calculations are made. The small line 
segments stemming from the circles indicate the changes 
of displacements at the locations. A scale of the 
displacement is given in the lower left corner and a 
scale of the map is given in the lower right corner. 

The thick line indicates the location of the earthquake 
fault. The boundaries of the descending slab and the 
low viscosity zone are also indicated on the plot, (b) 
The same plot at 3.6 years after the event. (c) The 
same plot at 11.0 years after the event. (d) The same 
plot as figure (c) with a reduced displacement scale. 
The boundaries of descending slab and low viscosity 
zones are also indicated in the figure. Notice that 
along the two sides of the dashed line, the 
displacements have opposite directions as if there were 
a dislocation there. 
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Figure 4.5.12 The viscosity ditribution for model 4M. The 
only difference from model 4L is that there is a longer 
descending slab. The zigzag contou** of the descending 
slab is due to the finite element discretization, but 
it has little effects on the vertical displacement on 
free surface. 

Figure 4.5.13 The vertical displacement change on free 

surface for 4M. The observed uplift data f>*om Bf'dwn et 
al are also superimposed on the figure. The model 
results are at 1.1, 3.6 and 11 years after the ev^nt. 
The observed uplifts are from 1965, 1968 and 1975. 

Figure 4.5.14 The vertical displacement change on free 

surface for layered model 4N. The observed uplift data 
from Brown et al are also superimposed on the figure. 
The model results are at 1.1, 3.6 and 11 years after 
the event. The observed uplifts a^e form 1965, 1968 
and 1975. 

Figure 4.6.1 The predicted vertical displacement change 
from Whittier to Anchorage f’-om model 4L . The curves 
indicated by 1, 2, 3, 4, 5 and 6 are at 1.1, 3.6, 11, 
13.2, 19.1 and 27.5 years after the event 
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Chapter 5 

Three dimensional analysis of dip slip earthquakes 
5 . 1 Introduction 

Great earthquakes may cause static crustal deformations 
detectable by modern instruments at teleseismic distances 
(Press 1965). This offers the possibility of studying 
earthquake mechanisms with teleseismic residual 
displacements, strains, and tilts. Great advances have been 
seen in new techniques of geodetic measurements in the last 
decade. In particular, lunar laser ranging (Bender and 
Silverberg 1975), satellite ranging (Agreen and Smith 1974, 
Smith et al 1979) and very long baseline interferometr y 
(Counselman 1973, Coates et al 1975, Niell et al 1979) all 
promise to provide distance measurements with the precision 
of a few cm per thousand km. A network of stations for zero 
frequency seismology will be feasible in the near future. 
Geodetic measurements can be used not only to study 
earthquake mechanisms but also to measure plate velocities 
in plate tectonics. The contributions to the fundamental 
understanding of geo^dynamics from such measurements will be 
tremendous. On the other hand, the modelling and 
interpretation of these data covering a great scope of space 
and time will certainly be a challenge for theoreticians. 

In anticipation of the coming data, in this, chapter 
we will make the first step in the modelling of the large 


scale time dependent problem with our finite element scheme. 
We will look into the problems of time dependent crustal 
movements at far away distances from the source and the 
stress diffusion due to a great dip slip event. 

The static displacements due to the 1964 Alaska 
earthquake were detected in Hawaii (Press 1965). The model 
of the I960 Chilean earthquake also indicated that the 
static effects caused by the earthquake were measurable in 
the plate interior (Richardson 1978b). However, the time 
dependent displacements due to great earthquakes at 
teleseismic distances have not been adequately studied. 

This information is essential in the deduction of plate 
velocities from geodetic data, since great plate boundary 
earthquakes may produce significant effects at the. stations 
of geodetic measurement. On the other hand, seismic 
activities in different regions are often correlated. 
Migration of earthquakes along the transform faults have 
been reported, as discussed in chapter 3. Migration and 
correlation of earthquakes in subduction zones have also 
been found (Fedotov 1965, Mogi 1968, 1969ab, Kelleher 1970, 
Sykes 1971, Kanamori 1972, Kelleher 1972, Mogi 1973, Utsu 
1974,1975, Shimazaki 1976,1978, Yamashina 1979). 

Correlation of global seismicity has been discovered using 
time series analysis (Chinnery and Landers 1975). 

Discussion of the migration of seismicity along subduction- 
zones needs a three dimensional model, because a two 
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dimensional plane strain model assumes that the fault length 
is infinite. In this chapter, we will investigate the 
possible triggering effects of a dip slip event. 

The modelling of large scale three dimensional problems 
is often hampered by the cost of computation. The number of 
operations in solving a system of linear equations in finite 
element analysis is proportional to NB , where N is the 
degree of freedom and B is the bandwidth of the stiffness 
matrix, which also often correlates with the size of the 
problem. When we discussed vertical strike slip earthquakes 
in chapter 3, the symmetry properties of the vertical strike 
slip fault were fully exploited. Only a quarter of the 
physical region was used in the computation. This 
constituted a great savings in computation cost. A dip slip 
event, except for the special case of a vertically dipping 
fault, can at best have a mirror symmetry. At least half of 
the physical region must be used in the computation. At the 
present time, there are no suitable geodetic data for 
comparison with the model results; we will use a coarse grid 
model to look at the general features of relaxation. 

However, the same scheme can be used with a finer grid model 
when data become available. 
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5.2 Model description 

We model the earthquake as a sudden slip (step 

function in time) on a fault surface. The medium is assumed 

to be elastic in bulk and viscoelastic in distortion. The 

12 2 

elastic properties are Young's modulus 2 x 10 dyne/ cm and 

Poisson's ratio 0,25. For the disturbances caused by the 

earthquake at far distances, the detailed structure in the 

source region is not important, and a layered model should 

be adequate. The model lithosphere has a thickness of 80 km 

25 

and a viscosity of 10 poise. From 80 km to 180 km the 

20 

viscosity is 10 poise, and below 180 km the viscosity is 
22 

10 poise. The schematic diagram of the fault and the 
viscosity model is given in figure 5.2.1. The fault is 500 
km long and 80 km deep, with dip angle 30 degrees. The fault 
slip is 10 m from 0 to 40 km deep, and it tapers off 
linearly from 10 m at 40 km deep to zero at 80 km deep; the 
fault slip also tapers off linearly to zero along the 50 km 
segments near the two fault edges along the strike 
direction. The finite element model uses half of the 
physical region. The top view of the finite element grid is 
given in figure 5*2.2. The plane y = 0 is the symmetry 
plane bisecting the fault. The fault intersects the free 
surface at coordinates x = 0 km and y = 0 to 250 km in the 
finite element model, as indicated by the dark line in the 
figure. The y displacement on the symmetry plane (y=0) is 
set to zero because of the symmetry. The lower boundary is 
assumed to be rigid while the side boundaries are assumed to 
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have spring boundary condition. The region is large enough 
that the external boundaries far away from the source have 
insignificant effects on the results. The physical model 
region is 3200 km long, 2200 km wide, and 700 km deep. Six 
grid surfaces and five layers of solid elements are used. 

The elements are eight node isoparametr ic hexahedrons with 
eight integration stations. There are 850 elements and 3^56 
degrees of freedom in the model. Admittedly this grid is 
too coarse to delineate the rapid variations near the 
source, however, it should be adequate to assess the general 
behavior away from the fault. 
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5.3 Model results 

We will only present the model results relevent to 
plate motion and stress diffusion. Figure 5.3.1 shows the 
surface horizontal displacement along the x axis (on the 
symmetry plane) at time = 0 and at 25 years after the 
event. More detailed time dependence can be seen from the 
plots of horizontal displacements vs. time at selected 
nodes. Figure 5.3.2 shows the surface horizontal 
displacements on the symmetry plane (y=0) at x = 0, -100, 
-400, -760, and -1000 km on the overiding plate side, and 
figure 5.3.3 shows the same quantities at x= 0, 100, 400, 
760, and 1000 km on the underthrusting plate side. The 
magnitude of the displacements gradually increase with time 
except, near the source. The increase is about 70 cm in four 
decades at 400 km away from the source, and it is less than 
20 cm at 1000 km away. The horizontal plate velocity due to 
the event averaged over four decades is less than 2 cm/year 
everywhere. However, the plate motion is nonuniform. Most 
of the time dependent adjustments take place within the 
first decade after the event. This can be seen more clearly 
in the plot of horizontal plate velocity in the x direction 
(V ) vs. time. Figure 5.3.4 gives V vs. time at the 

A A 

selected locations on the overiding plate side and figure 
5.3.5 gives on the und er thrusting plate side. The 
magnitude of V x decreases monotonically with time. At 
intermediate distance (x=400 km), V x is about 4 cm/year for 
the first several years after the event. This velocity is 
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comparable to the long term plate velocities in plate 
tectonic theory, and certainly is a measurable quantity in 
high precision geodetic measurements. On the other hand, V x 
is always less than 1 cm/year at 1000 km away from the 
source; this perturbation will be barely detectable with 
present day capabilities. 

Since the fault length is finite, the plate velocity 
caused by the earthquake varies with distance along the 
direction perpend icular to the symmetry plane. Figure 5.3.6 
shows V x along the line x = -400 km at y = 0, 120, and 250 
km, and figure 5.3*7 shows V x along the line x=-1000 km at y 
= 0, 210, and 540 km. V x in general decreases away from the 
symmetry plane (y=0). The edge effect of the source is 
significant at intermediate distance; along the line x = 

-400 km, the maximum V x decreases by 45 % at half fault 
length (250 km) away from the symmetry plane. This effect 
should be measurable but can not be described in two 
dimensional models. Along the line x = -1000 km, V x is 
again always less than 1 cm/year. 

Next we examine the stress diffusion caused by the 
earthquake. The largest stress component near the x axis is 
the horizontal normal stress along the x direction (S ). 
Figure 5.3.8 shows the stress component S at 20 km depth 

X X 

near the x axis at time = 0 and 25 years after the event. 
Positive S xx indicates horizontal tension and negative S xx 


* 


239 

indicates horizontal compression. The magnitude of S xx 
decreases with time near the source but increases with time 
at far distances, indicating a stress diffusion effect. The 
stress diffusion can be seen more clearly in the plot of 
stress vs. time. Figure 5.3.9 gives the stress component 
S at x = -60, -165, -385, -695, and -1175 km on the 

A A 

overiding plate side, and figure 5.3.10 gives S at x = 10, 

A A 

95, 225, .445, and 845 km on the underthrusting plate side. 
Near the source, the maxima of S occur at time = 0. At far 

A A 

away distances, the peaks of the stress arrive later, as if 
the stress is diffused away from the source region. If the 
source region is under horizontal compression before the 
earthquake, as expected for a thrust event, the thrust 
earthquake relieves the prestress near the source along the 
x axis. At some distance away from the source, the 
prestress state may be different from that of the source 
region. For example, the bending of the oceanic plate may 
produce tensile stress seawards of the thrust fault zone. A 
thrust fault event may reenforce the prestress there and 
trigger a normal fault event. 

The magnitude of the stress decreases rapidly with 
distance. At x= 1000 km, the magnitude of the stress is 
less than 2 bars. This stress probably can not produce 
significant physical effects by itself. However, the stress 
is long lasting and the curamulative effects due to many . 
events may be significant. 
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Earthquake migration along subduction zones has been 
reported (Mogi 1968, Kelleher 1970,1972). This may be 
related to the time dependent stress along the strike 
direction of a thrust earthquake. Figures 5.3.11 shows S xx 
vs. time at selected locations along the fault near the y 
axis. Figure 5.3.12 shows S xx at locations beyond the fault 
edge near the y axis. S xx changes sign near the fault edge 
along the strike direction. Along the fault zone S xx is 
positive (tension) while beyond the fault edge S xx is 
negative (compression). If the prestress near the source is 
horizontal compression, the thrust event relieves the 
prestress along the fault zone but reenforces the prestress 
in the zone beyond the fault edges along the strike 
direction. This is very similar to the stress diffusion 
associated with strike slip faults (chapter 3). Physically, 
this is not difficult to understand: if we consider the 
overiding plate to be stationary, the event tends to pull 
the underthrusting plate with it; this causes the 
compression in the adjacent regions. A great thrust 
earthquake thus increases the possibility of thrust 
earthquake occurrence in the zone adjacent to the source 
region along its strike direction. Notice that the 
magnitude of the compr essional stress due to the event 
increases with time; thus it is possible that the triggered 
event happens years after the triggering event. This is 
consistent with the earthquake migration phenomenon. 
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We would like to remind the reader that this model 
is a layered model. The effects due to the lateral 
viscosity variation in a subduction zone have not been 
included in the model. Those effects need to be calculated 
if detailed near source phenomena are discussed. 
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5.4 Discussion 

In this chapter, we made the first attempt to model 
the time dependent displacements and stress changes 
following a three dimensional dip slip earthquake. The main 
purpose of this study is to assess the time dependent 
effects of the earthquake on plate motion and stress 
diffusion. We used a simplified viscosity model and a 
coarse finite element grid. The effort needed in preparing 
a realistic three dimensional time dependent model of a 
subduction zone, including irregularly shaped descending 
slab and low viscosity zone, is huge; it should be reserved 
to the time when data warrant more detailed studies. 

However, the same computational scheme and programs can be 
used with little modification. 

The results in this chapter indicate that th& 
perturbation on the plate motion by a great thrust event is 
highly non-uniform in both space and time. For an 
earthquake with 10 m maximum slip, 500 km fault length, and 
30 degrees dip, on the symmetry plane the maximum horizontal 
surface velocity is about 4 cm/year at intermediate distance 
(about 400 km) from the source. This is comparable to long 
term average plate velocity. The maximum perturbation of 
plate velocity decreases to less than 1 cm/year at 1000 km 
away from the source. The results on the symmetry plane are 
similar to those of a two dimensional model. Figure 5.4.1 
gives the plate velocity at 380 and 900 km away from the 
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source for a two dimensional model. The model uses the 
finite element grid of figure 4.4.3 in chapter 4. The fault 
dip is 45 degrees and the slip is 10 m from surface to 40 km 
and tapers to zero at 80 km depth. The viscosity model is 
similar to that in this chapter. The maximum plate velocity 
is about 4 cm/ year at 380 km and 1 cm/year at 900 km. The 
maximum adjustments occur immediately after the event. 

These are similar to the three dimensional model. However, 
the decay of velocity with time is slower in the two 
dimensional model. Two dimensional model also cannot 
describe the effects off the symmetry plane. As we see in 
figure 5.3.6 and 5.3.7, the horizontal plate velocity varies 
considerably along the strike direction. 

The transient results in this chapter do not agree 
with the analyses using an elastic plate overlying a purely 
viscous fluid model (Elsasser 1967, Bott and Dean 1973, 
Anderson 1975, Melosh 1976). The difference is demonstrated 
in the following example. Figure 5.4.2 shows the schematic 
diagram for such a model: An elastic plate of thickness h^ 
simulates the lithosphere while a layer of viscous fluid of 
thickness simulates the asthenosphere . A step function 
displacement on the edge of the elastic plate simulates a 
plate boundary earthquake. If we neglect the bending in the 
lithospheric plate, then the solution for the horizontal 
displacement as a function of time and space is a 
complementary error function (Elsasser 1969). Assuming the 
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parameter values h^=80 km, h^^OO km, Young's modulus 

12 2 20 
2.0x10 dyne/cm , viscosity 10 poise, and initial 

displacement on the edge 10 m, the resulting plate 

velocities at 400 and 1000 km away are given in figure 

5.4.3. The transient velocity for the plate over viscous 

fluid model has a strong propagation effect. There is no 

instantaneous plate motion at any points other than at the 

source. The plate velocity reaches a maximum at a later 

time, then it decays to zero monotonically . This 

propagation effect of plate velocity does not exist in the 

viscoelastic models. The maximum adjustments in 

viscoelastic models happen right after the earthquake. The 

plate velocities are also much larger in the viscous fluid 

model than those in the viscoelastic model. There are two 

reasons for these differences: Viscous fluid cannot be 

deformed instantaneously by shearing, thus the plate 

velocity must be zero immediately after the event; the 

maximum plate velocity must occur later. In contrast, the 

viscoelastic models have instantaneous respone, and the 

maximum adjustments occur right after the event. This 

effect has also been pointed out by Savage and Prescott 

(1978a). The other reason is that dislocations are used to 

simulate the earthquakes in our viscoelastic models, while a 

step in displacement simulates the earthquake in the viscous 

fluid model. In general the effect of a dislocation decays 

with distance faster than a simple push or pull. This also 
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makes the plate velocities in the viscous fluid model 
larger. In the discussion of transient phenomena associated 
with an earthquake, the models assuming viscoelastic media 
and dislocation faults are more realistic than the simple 
viscous fluid model. 

The thrust event causes horizontal tension along the 
fault but horizontal compression in the regions beyond the 
fault edges along the strike direction. The tensile stress 
may trigger the normal fault events seawards of the thrust 
region; the compressional stress may trigger the thrust 
fault events in the region along the strike direction. The 
migration of thrust fault events along subduction zones has 
been observed in Aleutian, Kurile Islands, Japan, and South 
America (Fedotov 1965, Mogi 1968, Kezlleher 1970, 1972). The 
magnitude of the compressional stress caused by the thrust 
fault earthquake increases with time; thus the triggered 
events may happen years after the triggering event. This is 
consistent with the observed migration of earthquakes in 
subduction zones. 

Normal fault earthquakes triggered by thrust fault 
events have also been observed in subduction zones. A large 
normal fault event ( 111 ^= 7 ) occurred on March 30, 1965 beneath 
the Aleutian trench near Rat Island (Stauder 1968, Abe 
1972). This event was considered to be triggered by the 
thrust fault event on February 4, 1965 located at about 100 
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km away towards the island (Abe 1972). Normal fault events 
were also triggered in the 1938 Shioya-oki swarm of large 
earthquakes in the southern part of Japan trench (Abe 
1977b). In the Shioya-oki region, there had been no other 
earthquakes larger than magnitude 7.5 in historical time. A 
large earthquake of magnitude 7.4 occurred in May 1938; 
about 6 months later four earthquakes with magnitudes larger 
than 7 happened within 2 days. The fault mechanisms for 
these earthquakes were determined by Abe (1977) using 
seismic, tsunami and geodetic data. The first three events 
had shallow dipping thrust faults while the latter two had 
normal faults. The two normal fault events were located 
seawards of the three thrust fault events. The locations and 
occurence times suggested that the two normal fault events 
were triggered by the first three thrust fault events. 

Notice that near the source region, the tensile stress 
caused by a thrust fault event decays with time (figures 
5.3*9 and 5.3.10); thus the triggered normal faulting event 
should happen soon after the thrust fault event. In 
contrast, the migration of thrust fault events along the 
strike direction may happen years after the triggering 
event, because the triggering stress increases with time in 
this case (figure 5.3.12). 

The magnitude of the stress decreases rapidly with 
distance from the source. At 1000 km away, the magnitude of 
the stress is less than 2 bars. However, the stress at this 
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distance increases with time. For one event the physical 
effect is not significant. But the stress is long lasting 
and the cummulative effects due to .many events may be 
significant . 
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Figure Captions - Chapter 5 

Figure 5.2.1 The viscosity and fault of the three 

dimensional dip slip earthquake model in this chapter, 
(a) Cross sectional view of the plane (y=0) bisecting 
the fault. The fault has offset 10 m along the solid 
line segment and tapers off to zero at 80 km depth 
along the broken line segment. The fault dip is 30 
degrees. The viscosity has a layered structure. The 
numbers with exponents are the viscosity values in 
poise for each region, (b) Projection of the fault on 
the the surface. (z=0 plane). The fault intersects the 
surface at x=0 and yr-250 to 250 km. The fault has 
offset 10 m in the double hatched area and the offset 
tapers linearly in the single hatched area. 

Figure 5.2.2 Top view of the finite element grid used in 
the three dimensional model. The finite element model 
uses half of the physical model. The plane y=0 is the 
symmetry plane bisecting the physical region. The 
grids for surfaces at lower depths are similar to this 
top surface except for some shift to accomodate the 
fault configuration. 

Figure 5.3.1 The surface horizontal displacement along the 
symmetry plane (y=0). Curves 1 and 2 are for times 0 
and 25 years after the event. Positive values are in 
the positive X direction, i.e., toward ocean. 
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Figure 5.3*2 The surface horizontal displacement along the 
symmetry plane (y=0) vs. time at selected locations on 
the overiding plate side. Curves 1, 2, 3, 4, and 5 are 
for horizontal distances x = 0, -100, -400, -760, and 
-1000 km. Positive values are in positive X direction. 

Figure 5.3.3 The surface horizontal displacement along the 
symmetry plane (y=0) vs. time at selected locations on 
the underthrusting plate side. Curves 1, 2, 3, 4, and 
5 are for horizontal distances x= 0, 100, 400, 760, and 
1000 km. 

Figure 5.3*4 The surface horizontal velocity V x along the 
symmetry plane (y=0) vs. time at selected locations on 
the overiding plate side. Curves 1, 2, 3, 4, and 5 are 
for horizontal distances x= 0, -100, -400, -760, and 
-1000 km. 

Figure 5.3*5 The surface horizontal velocity V along the 
symmetry plane (y=0) vs. time at selected locations on 
the underthrusting plate side. Curves 1, 2, 3, 4, and 
5 are for horizontal distances x= 0, 100, 400, 760, and 
1000 km. 

Figure 5.3*6 The surface horizontal velocity V x along the 
line x = -400 km. The- curves 1, 2, and 3 are at y= 0, 

120, and 250 km. 
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Figure 5.3.7 The surface horizontal velocity V along the 

A 

line x=-1000 km. The curves 1, 2, and 3 are at y= 0, 
210, and 560 km. 

Figure 5.3.8 The normal stress perpendicular to the fault 
strike (S xx ) at 20 km depth near the x axis. Curves 1 
and 2 are for time = 0 and 25 years respectively. 
Positive stress indicates tension and negative stress 
indicates compression. 

Figure 5.3.9 S xx vs. time at selected locations near x 

axis at 20 km deep on the overiding plate side. Curves 
1, 2, 3, 4, and 5 are for x = -60, -165, -385, -695, 
and -1175 km. 

Figure 5.3.10 S xx vs. time at selected locations near x 
axis at 20 km deep on the underthrusting plate side. 
Curves 1, 2, 3, 4, and 5 are at x = 10, 95, 225, 445 
and 845 km. 

Figure 5.3.11 S xx vs. time at selected locations along the 
fault near the y axis (strike direction) at 20 km deep. 
The curves 1, 2, and 3 are at y = 25, 75, and 125 km. 

Figure 5.3*12 S xx vs. time at selected locations beyond the 
fault edge near the y axis (strike direction) at 20 km 
deep. The curves 4, 5, 6, and 7 are at y = 285, 380, 
510, and 670 km. 
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Figure 5.4.1 Typical time dependent behavior of plate 

velocity of a two dimensional viscoelastic model. The 

model uses the finite element grid of figure 4.4.3 in 

chapter 4. The fault dip is 45 degrees. The fault 

slip is 10 m from surface to 40 km deep; it tapers off 

to zero at 80 km deep. The viscosity structure is 

similar to figure 5.3.1 except that the layer with 
20 

viscosity 10 poise is 80 km thick. Curves 1 and 2 
are for 380 and 900 km away from the source. 

Figure 5.4.2 The schematic diagram of the plate over 

viscous fluid model. An elastic plate simulates the 
lithopshere. A viscous fluid layer simulates the 
asthenosphere . A sudden horizontal displacement on the 
edge of the elastic plat< ! ‘ simulates an earthquake. 

Figure 5.4.3 The horizontal plate velocities for the model 

in figure 5.4.2. The model parameters are: Young's 

12 2 

modulus in the lithosphere 2x10 dyne/cm , viscosity 

20 

in the asthenosphere 10 poise, h^ = 80 km, h 2 = 200 
km, and the sudden displacement on the edge of the 
lithosphere 10 m. Curves 1 and 2 are for x = 400 km 


and 1000 km. 
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Chapter 6 
Conclusions 

This thesis is mainly concerned with the problem of 
time dependent deformation and stress diffusion associated 
with earthquakes. A time dependent finite element scheme 
incorporating frontal solution and time stepping procedure 
was implemented. The efficiency of the scheme allowed us to 
model large scale problems on a moderate size computer. We 
applied the finite element analysis to models of two and 
three dimensional thrust earthquakes and three dimensional 
strike slip earthquakes. The finite element procedure 
permited us to model the viscosity distribution 
realistically. The modelling indicated that variations in 
viscosity distribution often produce different time 
dependent deformation patterns. The distinctive patterns of 
the time dependent deformation associated with an earthquake 
can be used to infer the viscosity structure of the fault 
zone. Stress relaxation following an earthquake is probably 
related to the seismicity in the adjacent regions. The 
stress diffusion effects in the models are consistent with 
the migration and correlation of seismicity in both 
transform fault zones and subduction zones. 

Models of three dimensional strike slip earthquakes 
indicated that lateral heterogeneities in viscosity across 
the fault zone can produce profound effects on the vertical 
deformation. If a low viscosity zone exists beneath the 
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fault zone, the increased viscosity away from the fault 
deters the material from moving away from the fault zone 
after the earthquake. This complicates the surface 
deformation. For example, the sudsiding region near fault 
for a layered model can become an uplift region with the 
laterally heterogeneous model. Thus the vertical 
deformation after a strike slip event, although smaller in 
magnitude than the horizontal deformation, is the most 
diagnostic information for the viscosity structure in a 
transform fault zone. The shear stress in the adjacent 
region along the strike direction increases with time due to 
the event. Therefore the possibility of earthquakes in the 
adjacent regions is enhanced in the years following the 
event. This stress diffusion may have caused the observed 
earthquake migration along the North Anatolian fault zone. 

For the postseismic vertical displacement following 
thrust events in simple layered models, the area near the 
lower fault tip subsides continuously while the adjacent 
areas are uplifted. Increasing the viscosity with depth 
reduces the subsidence and increases the uplifts. The 
relaxation in the simple layered model cannot explain the 
postseismic uplift data following the 196-4 Alaskan 
earthquake. Rather, certain models including a descending 
slab and a low viscosity wedge above the descending slab can 
reproduce the observed phenomenon. The descending slab acts 
as a guide for the induced viscous movement; it prevents the 
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material in the low viscosity wedge from going down and 
guides the movement sideways. The resulting deformation 
pattern in the downdip extension region of the earthquake 
fault is very similar to that of an aseismic slip, and the 
uplift data can be fit very well. These models, although 
still simple compared with the real world, probably contain 
the essential elements in a subduction zone. We believe 
that it is important to include these lateral variations of 
viscosity for modelling a subduction zone. 

In three dimensional model of a thrust earthquake, it 
was found that a great earthquake could produce measurable 
disturbance of plate velocity at intermediate distances. 
Beyond 1000 km away, the perturbation of plate velocity 
becomes smaller than 1 cm/year. The transient plate 
velocity of a viscoelastic model differs substantially from 
that of a simple elastic plate over viscous fluid models. 

A thrust fault earthquake causes horizontal compression 
in the adjacent regions along its strike direction. 

Therefore the possibility of another thrust earthquake in 
the adjacent region increases after the original event. 

This is consistent with the observation of the migration of 
great thrust events along subduction zones. The event 
causes horizontal tension in the adjacent regions along its 
dip direction; this is consistent with the observation of 
the normal fault events seawards of a thrust event. 
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With the available time dependent data today, often 
an observed phenomenon cannot be interpreted uniquely. 
However, the information on the time dependent deformation 
of the earth is being added at a rapid rate. There is no 
question that the study of time dependent behavior of the 
earth will provide significant new information on viscosity 
structure and earth dynamics. 
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